[卫生统计] 四参数拟合方程 的参数求解

[LanPst]

卫生统计学中 有对检验数据进行分析之后通过计算得到 样品的测量分析结果.
其中包含 四参数拟合方程 的方法(4 Parameter Logistic)找了很久都未找到.
不知道有没有人知道如何求解.


方程如下:(四参数为A B C D)
Y = (A - D) / (1 + (X / C) ^ B)) + D

给出数据集:
x y
0.00 1.4735
0.025 1.214
0.05 1.0143
0.1 0.7362
0.25 0.3907
0.75 0.2135

求最合适的四参数ABCD值是多少?

答案是 a = 1.47 b = 1.229 c = 0.085 d = 0.12

[LanPst]

自己顶 ~

[zhoushunda960]

Y = (A - D) / (1 + (X / C) ^ B)) + D 像这种方程，先归纳，然后改成y=ax^3 + bx^2 + cx + d的形式。
下面给你个算法，很好用，我一直用它做任意参数的拟合运算。

#define ABS(x) (x)>0?(x):-(x)
#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}

//求解线性方程组 A*x=B (n为矩阵维数)
void LinearSolve(double **A,double *B,double *x,int n)
{
	int i,j,k,ik;
	double mik,temp;

	double **a,*b;
	a=new double*[n];
	for(i=0;i<n;i++) a[i]=new double[n];

	b=new double[n];

	for(i=0;i<n;i++)
	{
		b[i]=B[i];
		for(j=0;j<n;j++) a[i][j]=A[i][j];
	}

	for(k=0;k<n;k++)
	{
		mik=-1;
		for(i=k;i<n;i++)
			if(ABS(a[i][k])>mik)
			{
				mik=ABS(a[i][k]);
				ik=i;
			}

		for(j=k;j<n;j++)
			SWAP(a[ik][j],a[k][j]);

		SWAP(b[k],b[ik]);
		b[k]/=a[k][k];
		for(i=n-1;i>=k;i--)
			a[k][i]/=a[k][k];

		for(i=k+1;i<n;i++)
		{
			b[i]-=a[i][k]*b[k];
			for(j=n-1;j>=k;j--)
				a[i][j]-=a[i][k]*a[k][j];
		}
	}
	for(i=n-1;i>=0;i--)
	{
		x[i]=b[i];
		for(j=i+1;j<n;j++)
			x[i]-=a[i][j]*x[j];
	}

	//删除内存
	delete b;
	for(i=0;i<n;i++) delete a[i];
}

//一般线性方程组拟合。n为参数个数,p为点数。
void linear(double **x,double *y,double *beta,int n,int p)
{
	double **a,*b;
	int i,j,k;

	//a:分配内存
	a=new double*[p];
	for(i=0;i<p;i++)
		a[i]=new double[p];

	//a:p*p矩阵
	for(i=0;i<p;i++)
		for(j=0;j<p;j++)
		{
			a[i][j]=0;
			for(k=0;k<n;k++)
				a[i][j]+=x[k][i]*x[k][j];	//计算内积
		}

	//b:1*p矩阵
	b=new double[p];
	for(i=0;i<p;i++)
	{
		b[i]=0;
		for(j=0;j<n;j++)
			b[i]+=x[j][i]*y[j];
	}

	LinearSolve(a,b,beta,p);	//求解beta

	//a:删除内存
	for(i=0;i<p;i++) delete a[i];

	delete b;
}

//多项式拟合：输入点坐标(x,y),输出多项式参数beta
//n为点数,p为多项式系数个数
void polyfit(double *x,double *y,double *beta,int n,int p)
{
	int i,j;
	double **xx=new double*[n];

	for(i=0;i<n;i++)
	{
		//多次拟合
		xx[i]=new double[p];
		for(j=0;j<p;j++)
		{
		    if(0==j) 
				xx[i][j]=1;
            else 
				xx[i][j]=pow(x[i],j);	//连续函数
		}
	}

	linear(xx,y,beta,n,p);

	for(i=0;i<n;i++) delete xx[i];

	delete xx;
}

不知道符合lz的要求。

[LanPst]

归纳是指什么？

B是X的指数，所以很难转化成明确的3，2，1次幂哦。。


[Quote=引用 4 楼 zhoushunda960 的回复:]
Y = (A - D) / (1 + (X / C) ^ B)) + D 像这种方程，先归纳，然后改成y=ax^3 + bx^2 + cx + d的形式。
下面给你个算法，很好用，我一直用它做任意参数的拟合运算。


C/C++ code

#define ABS(x) (x)>0?(x):-(x)
#define SWAP(a,b) {temp=(a);(a)=(b);(b)=te……
[/Quote]

[zh852]

属于线性方程 关于这个方面的问题 你可一使用lingo来解答 很方便的

[超级大笨狼]

这个东西还是Matlab的强项，秒杀是必须的。
建议去Matlab论坛问一下。

[LanPst]

额。。我现在要做一个拟和程序 - -
别的程序可以作验证，但计算过程估计是看不到了 - -

[Quote=引用 6 楼 zh852 的回复:]
属于线性方程 关于这个方面的问题 你可一使用lingo来解答 很方便的
[/Quote]

[LanPst]

- -我在MATLAB发贴 有20人查看却无一人回复
。压力。。。

[Quote=引用 7 楼 superdullwolf 的回复:]
这个东西还是Matlab的强项，秒杀是必须的。
建议去Matlab论坛问一下。
[/Quote]

[zhoushunda960]

上面的程序已经写得很清楚了。你直接调用polyfit函数就可以得出你要的结果。
double x（5）={0.00,0.025,0.05,0.1,0.25,0.75};
double y{5}={1.4735,1.214,1.0143,0.7362,0.3907,0.2135};
double beta(3);
polyfit(*x,*y,*beta,6,4);

最后beta里就是得到的A、B、C、D四个系数。

[zhoushunda960]

函数参数中的6是指数据的样品数据集元素个数，4只的是要得到的系数个数。

[zhoushunda960]

建议楼主稍微看一下数值分析，里面写得很详细。

[LanPst]

- -没人知道吗？

[yaoweijq]

找找matlab的资料吧
看能否秒杀

[LanPst]

我用别的软件验证够了，这个只能拟合多项式的方程，不能拟合幂函数的方程 - -
不过还是谢谢你。。。

[Quote=引用 10 楼 zhoushunda960 的回复:]
上面的程序已经写得很清楚了。你直接调用polyfit函数就可以得出你要的结果。
double x（5）={0.00,0.025,0.05,0.1,0.25,0.75};
double y{5}={1.4735,1.214,1.0143,0.7362,0.3907,0.2135};
double beta(3);
polyfit(*x,*y,*beta,6,4);

最后beta里就……
[/Quote]

[LanPst]

不知道你的程序算出来结果是多少，
我算的话只有第一个和答案一样，其他三个都不一样。
-10.2
29.6
-24.5
我把程序改成C#的了.不知道程序是不是有问题.你算的是这个答案吗?
另:这个方法是不是拟合成多项式的?
形如:ax^3+bx^2+cx+d ?
这个和要求的不一样哦.因为 要求中参数b是x的指数,而不是系数.

[Quote=引用 10 楼 zhoushunda960 的回复:]
上面的程序已经写得很清楚了。你直接调用polyfit函数就可以得出你要的结果。
double x（5）={0.00,0.025,0.05,0.1,0.25,0.75};
double y{5}={1.4735,1.214,1.0143,0.7362,0.3907,0.2135};
double beta(3);
polyfit(*x,*y,*beta,6,4);

最后beta里就……
[/Quote]

[cnmhx]

这个问题得用非线性拟合。
譬如买夸特/共轭梯度。。。等方法。

[chiyuq]

我用C++算法递归拟合结果
a0+(a1-a0)/(1+pow(x/a2,a3))
.04663,1.4646, .0922, 1.0135 基本一致,可能精度误差

[hyphen0104]

如何解出来的呢？

[Quote=引用 16 楼 chiyuq 的回复:]

我用C++算法递归拟合结果
a0+(a1-a0)/(1+pow(x/a2,a3))
.04663,1.4646, .0922, 1.0135 基本一致,可能精度误差
[/Quote]

[火星大能猫]

表示关注这个问题

[火星大能猫]

[Quote=引用 16 楼 的回复:]

我用C++算法递归拟合结果
a0+(a1-a0)/(1+pow(x/a2,a3))
.04663,1.4646, .0922, 1.0135 基本一致,可能精度误差
[/Quote]
是啊.怎么算的啊???/