2013-02-01 32 views
2

对于几天我一直对这个问题,我被困...曲线拟合R中的这些数据?

已经执行了数R中蒙特卡洛模拟的,其给出了每个输入x输出y和有明确x和y之间有一些简单的关系,所以我想识别公式及其参数。但我似乎无法为'低x'和'高x'两个系列获得良好的整体配合,例如使用对数是这样的:

dat = data.frame(x=x, y=y) 
fit = nls(y~a*log10(x)+b, data=dat, start=list(a=-0.8,b=-2), trace=TRUE) 

我也试图适应(LOG10(X),10^Y)代替,它给出了一个不错的选择,但相反转变不适合(X,Y),非常好。

任何人都可以解决这个问题吗?

请解释您是如何找到解决方案的。

谢谢!

编辑:

感谢所有的快速反馈!

我不知道我模拟的理论模型,所以我没有比较的基础。我根本不知道x和y之间的真实关系。顺便说一句,我不是统计学家。

基础模型是一种随机反馈增长模型。我的目标是在给定输入x> 0的情况下确定长期增长率g,因此系统的输出在每次迭代中按照1 + g的指数增长。系统在每次迭代时都会根据系统的规模进行随机生产,其中一部分生产是输出的,其余部分则保存在由另一个随机变量确定的系统中。从MC模拟中我发现系统输出的增长率对于我测试的每个x都是对数正态分布的,并且数据系列中的y是增长率g的对数。当x走向无穷时,g趋于零。当x走向零时,g走向无穷大。

我想要一个可以从x计算y的函数。实际上我只需要一个低x的函数,比如在0到10的范围内。我可以很好地拟合y = 1.556 * x^-0.4 -3.58,但是它对于大x很不适合。我想要一个对所有x> 0都是通用的函数。我也尝试过Spacedman的poly fit(谢谢!),但它在x = 1到6的关键范围内不够好。

任何想法?

编辑2:

我已经尝试格罗腾迪克多一些,也有详细的建议,经过一番考虑,我决定,因为我没有超过选择一个功能提供了理论依据(感谢!)另一个,我很可能只对1到6之间的x值感兴趣,所以我应该使用一个很好的简单函数。所以我只是使用了y〜a * x^b + c并且表示它不适合高x。当论文初稿完成后,我可能会再次寻求社区的帮助。一旦你看到蒙特卡罗模型,或许你们中的一个人可以发现x和y之间的理论关系。再次

谢谢!

低X系列:

 x   y 
1 0.2 -0.7031864 
2 0.3 -1.0533648 
3 0.4 -1.3019655 
4 0.5 -1.4919278 
5 0.6 -1.6369545 
6 0.7 -1.7477481 
7 0.8 -1.8497117 
8 0.9 -1.9300209 
9 1.0 -2.0036842 
10 1.1 -2.0659970 
11 1.2 -2.1224324 
12 1.3 -2.1693986 
13 1.4 -2.2162889 
14 1.5 -2.2548485 
15 1.6 -2.2953162 
16 1.7 -2.3249750 
17 1.8 -2.3570141 
18 1.9 -2.3872684 
19 2.0 -2.4133978 
20 2.1 -2.4359624 
21 2.2 -2.4597122 
22 2.3 -2.4818787 
23 2.4 -2.5019371 
24 2.5 -2.5173966 
25 2.6 -2.5378936 
26 2.7 -2.5549524 
27 2.8 -2.5677939 
28 2.9 -2.5865958 
29 3.0 -2.5952558 
30 3.1 -2.6120607 
31 3.2 -2.6216831 
32 3.3 -2.6370452 
33 3.4 -2.6474608 
34 3.5 -2.6576862 
35 3.6 -2.6655606 
36 3.7 -2.6763866 
37 3.8 -2.6881303 
38 3.9 -2.6932310 
39 4.0 -2.7073198 
40 4.1 -2.7165035 
41 4.2 -2.7204063 
42 4.3 -2.7278532 
43 4.4 -2.7321731 
44 4.5 -2.7444773 
45 4.6 -2.7490365 
46 4.7 -2.7554178 
47 4.8 -2.7611471 
48 4.9 -2.7719188 
49 5.0 -2.7739299 
50 5.1 -2.7807113 
51 5.2 -2.7870781 
52 5.3 -2.7950429 
53 5.4 -2.7975677 
54 5.5 -2.7990999 
55 5.6 -2.8095955 
56 5.7 -2.8142453 
57 5.8 -2.8162046 
58 5.9 -2.8240594 
59 6.0 -2.8272394 
60 6.1 -2.8338866 
61 6.2 -2.8382038 
62 6.3 -2.8401935 
63 6.4 -2.8444915 
64 6.5 -2.8448382 
65 6.6 -2.8512086 
66 6.7 -2.8550240 
67 6.8 -2.8592950 
68 6.9 -2.8622220 
69 7.0 -2.8660817 
70 7.1 -2.8710430 
71 7.2 -2.8736998 
72 7.3 -2.8764701 
73 7.4 -2.8818748 
74 7.5 -2.8832696 
75 7.6 -2.8833351 
76 7.7 -2.8891867 
77 7.8 -2.8926849 
78 7.9 -2.8944987 
79 8.0 -2.8996780 
80 8.1 -2.9011012 
81 8.2 -2.9053911 
82 8.3 -2.9063661 
83 8.4 -2.9092228 
84 8.5 -2.9135426 
85 8.6 -2.9101730 
86 8.7 -2.9186316 
87 8.8 -2.9199631 
88 8.9 -2.9199856 
89 9.0 -2.9239220 
90 9.1 -2.9240167 
91 9.2 -2.9284608 
92 9.3 -2.9294951 
93 9.4 -2.9310985 
94 9.5 -2.9352370 
95 9.6 -2.9403694 
96 9.7 -2.9395336 
97 9.8 -2.9404153 
98 9.9 -2.9437564 
99 10.0 -2.9452175 

(高)x系列:

   x   y 
1 2.000000e-01 -0.701301 
2 2.517851e-01 -0.907446 
3 3.169786e-01 -1.104863 
4 3.990525e-01 -1.304556 
5 5.023773e-01 -1.496033 
6 6.324555e-01 -1.674629 
7 7.962143e-01 -1.842118 
8 1.002374e+00 -1.998864 
9 1.261915e+00 -2.153993 
10 1.588656e+00 -2.287607 
11 2.000000e+00 -2.415137 
12 2.517851e+00 -2.522978 
13 3.169786e+00 -2.621386 
14 3.990525e+00 -2.701105 
15 5.023773e+00 -2.778751 
16 6.324555e+00 -2.841699 
17 7.962143e+00 -2.900664 
18 1.002374e+01 -2.947035 
19 1.261915e+01 -2.993301 
20 1.588656e+01 -3.033517 
21 2.000000e+01 -3.072003 
22 2.517851e+01 -3.102536 
23 3.169786e+01 -3.138539 
24 3.990525e+01 -3.167577 
25 5.023773e+01 -3.200739 
26 6.324555e+01 -3.233111 
27 7.962143e+01 -3.259738 
28 1.002374e+02 -3.291657 
29 1.261915e+02 -3.324449 
30 1.588656e+02 -3.349988 
31 2.000000e+02 -3.380031 
32 2.517851e+02 -3.405850 
33 3.169786e+02 -3.438225 
34 3.990525e+02 -3.467420 
35 5.023773e+02 -3.496026 
36 6.324555e+02 -3.531125 
37 7.962143e+02 -3.558215 
38 1.002374e+03 -3.587526 
39 1.261915e+03 -3.616800 
40 1.588656e+03 -3.648891 
41 2.000000e+03 -3.684342 
42 2.517851e+03 -3.716174 
43 3.169786e+03 -3.752631 
44 3.990525e+03 -3.786956 
45 5.023773e+03 -3.819529 
46 6.324555e+03 -3.857214 
47 7.962143e+03 -3.899199 
48 1.002374e+04 -3.937206 
49 1.261915e+04 -3.968795 
50 1.588656e+04 -4.015991 
51 2.000000e+04 -4.055811 
52 2.517851e+04 -4.098894 
53 3.169786e+04 -4.135608 
54 3.990525e+04 -4.190248 
55 5.023773e+04 -4.237104 
56 6.324555e+04 -4.286103 
57 7.962143e+04 -4.332090 
58 1.002374e+05 -4.392748 
59 1.261915e+05 -4.446233 
60 1.588656e+05 -4.497845 
61 2.000000e+05 -4.568541 
62 2.517851e+05 -4.628460 
63 3.169786e+05 -4.686546 
64 3.990525e+05 -4.759202 
65 5.023773e+05 -4.826938 
66 6.324555e+05 -4.912130 
67 7.962143e+05 -4.985855 
68 1.002374e+06 -5.070668 
69 1.261915e+06 -5.143341 
70 1.588656e+06 -5.261585 
71 2.000000e+06 -5.343636 
72 2.517851e+06 -5.447189 
73 3.169786e+06 -5.559962 
74 3.990525e+06 -5.683828 
75 5.023773e+06 -5.799319 
76 6.324555e+06 -5.929599 
77 7.962143e+06 -6.065907 
78 1.002374e+07 -6.200967 
79 1.261915e+07 -6.361633 
80 1.588656e+07 -6.509538 
81 2.000000e+07 -6.682960 
82 2.517851e+07 -6.887793 
83 3.169786e+07 -7.026138 
84 3.990525e+07 -7.227990 
85 5.023773e+07 -7.413960 
86 6.324555e+07 -7.620247 
87 7.962143e+07 -7.815754 
88 1.002374e+08 -8.020447 
89 1.261915e+08 -8.229911 
90 1.588656e+08 -8.447927 
91 2.000000e+08 -8.665613 
+0

您可以制作一个可重复使用的小例子吗? http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example –

+1

你能提供更多的信息吗?什么是x和y?非线性模型应该基于物理,生物,化学......考虑。顺便说一句,你显示的模型是线性的。不需要使用'nls'。 – Roland

+0

“低x”和“高x”是什么意思?它们的范围重叠。 –

回答

7

倒退X/Y与X绘制yx的低数据和逛了一下,似乎x/y被近似线性x打了这么尝试回归x/yx这给了我们仅仅基于关系两个参数:

y = x/(a + b * x) 

其中a和b是回归系数。

> lm(x/y ~ x, lo.data) 

Call: 
lm(formula = x/y ~ x, data = lo.data) 

Coefficients: 
(Intercept)   x 
    -0.1877  -0.3216 

MM.2上述可转化为在DRCř包MM.2模型。如下所示,该型号具有较高的R 。另外,我们计算了AIC,我们可以用它来比较其他车型(越低越好):

> library(drc) 
> fm.mm2 <- drm(y ~ x, data = lo.data, fct = MM.2()) 
> cor(fitted(fm.mm2), lo.data$y)^2 
[1] 0.9986303 
> AIC(fm.mm2) 
[1] -535.7969 

CRS.6这表明我们尝试一些其他的DRC模式和那些我们试图CRS。 6具有特别低的AIC,似乎合身视觉:

> fm.crs6 <- drm(y ~ x, data = lo.data, fct = CRS.6()) 
> AIC(fm.crs6) 
[1] -942.7866 
> plot(fm.crs6) # see output below 

这给我们带来了多款车型,我们可以从2参数MM.2模型,它是不是一个合适的好使用(根据AIC)作为CRS.6,但仍然非常适合,并具有只有两个参数或6参数CRS.6型号以其卓越的AIC。请注意,AIC已经对具有更多参数的模型进行惩罚,因此拥有更好的AIC不是具有更多参数的结果。

其他如果其认为,低和高应具有相同的标准格式,然后找到一个单一的模式形状配合低和高以及可作为另一个标准采摘模式的形式。除了drc模型外,在(2.1),(2.2),(2.3)和(2.4)Akbar et al, IRJFE, 2010中也有一些产量密度模型,它们看起来与可以尝试的MM.2模型相似。

screenshot

更新时间:整顿这个周围的DRC包。

+0

这是很久以前,我不记得的细节了。但我被其他人要求将其中的一个标记为答案,所以我标记了这一个,并将其中一个标记为另一个。再次感谢您的帮助。 – questiondude

9

没有底层程序的想法,你可能也只是适合与许多组件作为多项式你喜欢。你似乎没有测试一个假设(例如,引力强度是与距离相反的平方反比),所以你可以为功能表单捕鱼,数据不太可能告诉你哪一个是“正确的”。

所以,如果我看了你的资料与x和y分量数据帧我可以这样做:

data$lx=log(data$x) 
plot(data$lx,data$y) # needs at least a cubic polynomial 
m1 = lm(y~poly(lx,3),data=data) # fit a cubic 
points(data$lx,fitted(m1),pch=19) 

和拟合点非常接近。将多项式的等级从3更改为7,并且这些点是相同的。这是否意味着你的Y值真的来自X值的7次多项式?不可以。但是你有一条通过点的曲线。

在这个尺度上,你可能只是用一条直线连接相邻点,你的情节就非常平滑。但是,没有Y为什么依赖于X的基本理论(如反平方律或指数增长等),所有你正在做的就是加入点,并且有无限的方式来做到这一点。

+0

如果我找到另一个模型,例如在log(x)的函数中,例如使用相同的R square,我该如何选择2个模型? – agstudy

+3

@agstudy,除非您使用数据验证理论关系,否则您使用的是什么拟合函数并不重要。那就是找到合适的功能的原因。是预测不在原始数据集中的x值的输出。 (关于推断的常见可怕警告在这里)。 –

+0

@CarlWitthoft谢谢! – agstudy