I have following data
Rand_Data = data.frame(x = c( 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.,
11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21.,
22., 23., 24., 25., 26., 27., 28., 29., 30., 31., 32.,
33., 34., 35., 36., 37., 38., 39., 40., 41., 42., 43.,
44., 45., 46., 47., 48., 49., 50., 51., 52., 53., 54.,
55., 56., 57., 58., 59., 60., 61., 62., 63., 64., 65.,
66., 67., 68., 69., 70., 71., 72., 73., 74., 75., 76.,
77., 78., 79., 80., 81., 82., 83., 84., 85., 86., 87.,
88., 89., 90., 91., 92., 93., 94., 95., 96., 97., 98.,
99., 100., 101., 102., 103., 104., 105., 106., 107., 108., 109.,
110., 111., 112., 113., 114., 115., 116., 117., 118., 119., 120.,
121., 122., 123., 124., 125., 126., 127., 128., 129., 130., 131.,
132., 133., 134., 135., 136., 137., 138., 139., 140., 141., 142.,
143., 144., 145., 146., 147., 148., 149., 150., 151., 152., 153.,
154., 155., 156., 157., 158., 159., 160., 161., 162., 163., 164.,
165., 166., 167., 168., 169., 170., 171., 172., 173., 174., 175.,
176., 177., 178., 179., 180., 181., 182., 183., 184., 185., 186.,
187., 188., 189., 190., 191., 192., 193., 194., 195., 196., 197.,
198., 199), Val = c( 0.00000000e+00, 7.75383362e-02, -1.34275505e+00, -1.65497205e+00,
-1.18564171e+00, 1.77660878e+00, 1.84593088e+00, 1.71963694e+00,
8.21964340e-01, 1.74041137e+00, 2.64637861e+00, 2.06965750e+00,
1.02828336e+00, -5.00395684e-01, -6.93668766e-02, 4.22030196e-01,
1.16198422e-01, -5.11722287e-04, -6.57944961e-01, -3.74434142e-01,
-1.15827213e+00, -1.45830766e+00, -4.90458625e-01, 2.47042268e-02,
6.93327973e-01, 1.30938478e+00, -4.06142808e-01, -4.61168289e-01,
-6.72137769e-01, -4.84899922e-01, -8.29915365e-01, -2.41828942e+00,
-1.20765939e+00, 5.76353406e-01, 2.61955788e+00, 2.98795810e+00,
4.16904195e+00, 2.93858016e+00, 5.35527961e+00, 5.11840962e+00,
5.51236637e+00, 5.88671081e+00, 5.38415506e+00, 6.32405264e+00,
5.65838016e+00, 5.07960498e+00, 5.70407219e+00, 7.29978260e+00,
6.08456290e+00, 7.43051089e+00, 8.55807191e+00, 8.52981505e+00,
8.23329180e+00, 9.10450577e+00, 9.85288281e+00, 9.62989249e+00,
8.69069599e+00, 8.34052238e+00, 7.19647733e+00, 7.24917949e+00,
6.63422915e+00, 4.95911711e+00, 3.49654704e+00, 1.85980401e+00,
3.08957065e+00, 4.24267262e+00, 2.38508285e+00, 1.02089838e+00,
2.19670308e-01, 1.53858340e+00, 2.41653182e+00, 3.69861411e+00,
2.87607870e+00, 3.21819643e+00, 2.07359072e+00, 1.13000323e+00,
1.17047772e+00, 1.55913045e+00, 8.41514076e-01, 2.23618634e+00,
2.99905073e+00, 3.21259672e+00, 4.39274819e+00, 1.99515000e+00,
1.16175606e+00, 1.70758890e+00, 9.74160444e-01, -3.03140114e-02,
1.17659213e+00, 3.07232323e-01, 1.08034381e+00, 1.29199568e+00,
2.45607387e-01, -5.21021566e-01, -4.34639402e-01, -9.60127469e-01,
-1.03258482e+00, -1.43325172e+00, -2.28217823e+00, -2.09923712e+00,
-2.06051134e+00, -2.19767637e+00, -2.16617856e+00, -2.45811083e+00,
-2.79865309e+00, -3.70529459e+00, -2.77461502e+00, -3.27216972e+00,
-4.79844244e+00, -3.42557001e+00, -4.14314358e+00, -5.25738457e+00,
-3.49827107e+00, -3.55643468e+00, -5.45626633e+00, -5.33106305e+00,
-5.25332322e+00, -4.94845946e+00, -5.04795700e+00, -4.75058961e+00,
-3.37815392e+00, -3.24099974e+00, -4.26774262e+00, -4.35293178e+00,
-3.70409130e+00, -4.85062599e+00, -4.46039428e+00, -4.73719944e+00,
-5.53392228e+00, -6.15049133e+00, -6.10977286e+00, -4.87411131e+00,
-4.49530898e+00, -4.66802205e+00, -6.98105638e+00, -5.77006061e+00,
-6.27152149e+00, -6.07877879e+00, -5.09850276e+00, -5.68772276e+00,
-5.16446602e+00, -3.58115497e+00, -3.21209072e+00, -3.61179089e+00,
-4.40163561e+00, -3.14328137e+00, -3.80444247e+00, -3.91338929e+00,
-3.98395412e+00, -4.71086150e+00, -2.98308294e+00, -2.51608443e+00,
-2.86848586e+00, -3.27653288e+00, -3.17684302e+00, -2.88103208e+00,
-4.90737736e+00, -4.27541967e+00, -2.92691377e+00, -2.19142665e+00,
-3.08185258e+00, -2.54393704e+00, -2.42485142e+00, -1.90176762e+00,
-6.52372476e-01, -1.04169430e+00, -5.79717797e-01, -4.03696696e-01,
5.92495881e-01, -2.14694942e+00, -1.42679299e+00, -9.81665751e-01,
-3.51830293e+00, -2.24672044e+00, -2.23578918e+00, -9.44689326e-01,
-1.55913248e+00, -6.17618042e-01, -9.56082794e-01, -1.38326771e+00,
-1.55108670e+00, -5.74061887e-01, 1.89141995e-01, 8.11432327e-01,
1.99574115e+00, 1.72686176e+00, 7.17859558e-01, 5.03518805e-01,
1.91256417e+00, 9.59411484e-01, 5.28613488e-01, 1.17086649e+00,
1.88075943e+00, 2.41999514e+00, 2.51115680e+00, 4.77976432e+00,
6.14323607e+00, 6.42564314e+00, 4.82006996e+00, 4.63386067e+00))
And, below is my function (for simplicity, I just consider a Cubic spline fitting)
MyFN = function(a, b, c, d) return(a * Rand_Data$x^3 + b * Rand_Data$x^2 + c * Rand_Data$x + d)
Now I want to estimate the parameters a, b, c, d that minimises the sum of absolute difference between Rand_Data$Val and the fitted value with below constraints
MyFN should be less than 2 i.e. MyFN[100] < 2MyFN$Val should be less that 0.001I wonder how can I setup optim() function so that I can estimate a,b,c,d
Here, my data and function and just for example, so I am looking for some generic approach to optimally estimate the parameters.
You are solving constrained optimization problem, where nonlinear constraints are presented. I guess you could use fmincon from package pracma, e.g.,
library(pracma)
x <- Rand_Data$x
y <- Rand_Data$Val
eps <- 1e-4
X <- outer(x, 0:3, `^`)
# define functions
f <- function(theta) {
sum(abs(X %*% theta - y))
}
# define constraints
A <- X[100, , drop = FALSE]
b <- 2 - eps
hin <- function(theta) {
abs(sum(X %*% theta) - sum(y)) - 1e-3 + eps
}
# use fmincon
p <- fmincon(rep(0, 4), f, A = A, b = b, hin = hin)
sol <- p$par
res <- transform(
Rand_Data,
Val_est = X %*% sol
)
# verify the solution
with(
res,
abs(sum(Val) - sum(Val_est)) < 1e-3 & Val_est[100] < 2
)
and you will obtain the solution (d,c,b,a) like below
> sol
[1] -1.281027e+03 5.661665e+01 -6.601774e-01 2.214139e-03
and the estimated Val, namely, Val_est looks like
> res
x Val Val_est
1 0 0.0000000000 -1281.027073
2 1 0.0775383362 -1225.068384
3 2 -1.3427550500 -1170.416765
4 3 -1.6549720500 -1117.058931
5 4 -1.1856417100 -1064.981597
6 5 1.7766087800 -1014.171479
7 6 1.8459308800 -964.615291
8 7 1.7196369400 -916.299749
9 8 0.8219643400 -869.211568
10 9 1.7404113700 -823.337462
11 10 2.6463786100 -778.664148
12 11 2.0696575000 -735.178341
13 12 1.0282833600 -692.866755
14 13 -0.5003956840 -651.716105
15 14 -0.0693668766 -611.713108
16 15 0.4220301960 -572.844477
17 16 0.1161984220 -535.096929
18 17 -0.0005117223 -498.457178
19 18 -0.6579449610 -462.911940
20 19 -0.3744341420 -428.447929
21 20 -1.1582721300 -395.051861
22 21 -1.4583076600 -362.710452
23 22 -0.4904586250 -331.410415
24 23 0.0247042268 -301.138467
25 24 0.6933279730 -271.881322
26 25 1.3093847800 -243.625697
27 26 -0.4061428080 -216.358305
28 27 -0.4611682890 -190.065862
29 28 -0.6721377690 -164.735083
30 29 -0.4848999220 -140.352684
31 30 -0.8299153650 -116.905379
32 31 -2.4182894200 -94.379884
33 32 -1.2076593900 -72.762914
34 33 0.5763534060 -52.041184
35 34 2.6195578800 -32.201409
36 35 2.9879581000 -13.230305
37 36 4.1690419500 4.885414
38 37 2.9385801600 22.159033
39 38 5.3552796100 38.603835
40 39 5.1184096200 54.233106
41 40 5.5123663700 69.060131
42 41 5.8867108100 83.098195
43 42 5.3841550600 96.360582
44 43 6.3240526400 108.860577
45 44 5.6583801600 120.611465
46 45 5.0796049800 131.626531
47 46 5.7040721900 141.919061
48 47 7.2997826000 151.502337
49 48 6.0845629000 160.389646
50 49 7.4305108900 168.594273
51 50 8.5580719100 176.129501
52 51 8.5298150500 183.008616
53 52 8.2332918000 189.244904
54 53 9.1045057700 194.851647
55 54 9.8528828100 199.842133
56 55 9.6298924900 204.229644
57 56 8.6906959900 208.027466
58 57 8.3405223800 211.248885
59 58 7.1964773300 213.907184
60 59 7.2491794900 216.015649
61 60 6.6342291500 217.587564
62 61 4.9591171100 218.636214
63 62 3.4965470400 219.174885
64 63 1.8598040100 219.216860
65 64 3.0895706500 218.775426
66 65 4.2426726200 217.863865
67 66 2.3850828500 216.495465
68 67 1.0208983800 214.683508
69 68 0.2196703080 212.441281
70 69 1.5385834000 209.782067
71 70 2.4165318200 206.719152
72 71 3.6986141100 203.265821
73 72 2.8760787000 199.435358
74 73 3.2181964300 195.241049
75 74 2.0735907200 190.696177
76 75 1.1300032300 185.814029
77 76 1.1704777200 180.607888
78 77 1.5591304500 175.091040
79 78 0.8415140760 169.276769
80 79 2.2361863400 163.178360
81 80 2.9990507300 156.809098
82 81 3.2125967200 150.182269
83 82 4.3927481900 143.311156
84 83 1.9951500000 136.209045
85 84 1.1617560600 128.889220
86 85 1.7075889000 121.364966
87 86 0.9741604440 113.649568
88 87 -0.0303140114 105.756312
89 88 1.1765921300 97.698481
90 89 0.3072323230 89.489360
91 90 1.0803438100 81.142235
92 91 1.2919956800 72.670390
93 92 0.2456073870 64.087111
94 93 -0.5210215660 55.405681
95 94 -0.4346394020 46.639386
96 95 -0.9601274690 37.801510
97 96 -1.0325848200 28.905339
98 97 -1.4332517200 19.964157
99 98 -2.2821782300 10.991249
100 99 -2.0992371200 1.999900
101 100 -2.0605113400 -6.996605
102 101 -2.1976763700 -15.984982
103 102 -2.1661785600 -24.951945
104 103 -2.4581108300 -33.884210
105 104 -2.7986530900 -42.768492
106 105 -3.7052945900 -51.591507
107 106 -2.7746150200 -60.339968
108 107 -3.2721697200 -69.000592
109 108 -4.7984424400 -77.560094
110 109 -3.4255700100 -86.005188
111 110 -4.1431435800 -94.322590
112 111 -5.2573845700 -102.499015
113 112 -3.4982710700 -110.521179
114 113 -3.5564346800 -118.375796
115 114 -5.4562663300 -126.049582
116 115 -5.3310630500 -133.529251
117 116 -5.2533232200 -140.801520
118 117 -4.9484594600 -147.853103
119 118 -5.0479570000 -154.670714
120 119 -4.7505896100 -161.241071
121 120 -3.3781539200 -167.550887
122 121 -3.2409997400 -173.586878
123 122 -4.2677426200 -179.335758
124 123 -4.3529317800 -184.784244
125 124 -3.7040913000 -189.919050
126 125 -4.8506259900 -194.726892
127 126 -4.4603942800 -199.194484
128 127 -4.7371994400 -203.308542
129 128 -5.5339222800 -207.055781
130 129 -6.1504913300 -210.422916
131 130 -6.1097728600 -213.396663
132 131 -4.8741113100 -215.963736
133 132 -4.4953089800 -218.110850
134 133 -4.6680220500 -219.824722
135 134 -6.9810563800 -221.092065
136 135 -5.7700606100 -221.899595
137 136 -6.2715214900 -222.234028
138 137 -6.0787787900 -222.082078
139 138 -5.0985027600 -221.430461
140 139 -5.6877227600 -220.265892
141 140 -5.1644660200 -218.575086
142 141 -3.5811549700 -216.344757
143 142 -3.2120907200 -213.561623
144 143 -3.6117908900 -210.212396
145 144 -4.4016356100 -206.283793
146 145 -3.1432813700 -201.762529
147 146 -3.8044424700 -196.635319
148 147 -3.9133892900 -190.888879
149 148 -3.9839541200 -184.509922
150 149 -4.7108615000 -177.485165
151 150 -2.9830829400 -169.801323
152 151 -2.5160844300 -161.445110
153 152 -2.8684858600 -152.403242
154 153 -3.2765328800 -142.662435
155 154 -3.1768430200 -132.209403
156 155 -2.8810320800 -121.030861
157 156 -4.9073773600 -109.113525
158 157 -4.2754196700 -96.444110
159 158 -2.9269137700 -83.009331
160 159 -2.1914266500 -68.795903
161 160 -3.0818525800 -53.790541
162 161 -2.5439370400 -37.979961
163 162 -2.4248514200 -21.350878
164 163 -1.9017676200 -3.890006
165 164 -0.6523724760 14.415939
166 165 -1.0416943000 33.580241
167 166 -0.5797177970 53.616187
168 167 -0.4036966960 74.537059
169 168 0.5924958810 96.356145
170 169 -2.1469494200 119.086727
171 170 -1.4267929900 142.742091
172 171 -0.9816657510 167.335522
173 172 -3.5183029300 192.880305
174 173 -2.2467204400 219.389724
175 174 -2.2357891800 246.877065
176 175 -0.9446893260 275.355611
177 176 -1.5591324800 304.838649
178 177 -0.6176180420 335.339463
179 178 -0.9560827940 366.871337
180 179 -1.3832677100 399.447557
181 180 -1.5510867000 433.081407
182 181 -0.5740618870 467.786172
183 182 0.1891419950 503.575137
184 183 0.8114323270 540.461587
185 184 1.9957411500 578.458807
186 185 1.7268617600 617.580081
187 186 0.7178595580 657.838695
188 187 0.5035188050 699.247932
189 188 1.9125641700 741.821079
190 189 0.9594114840 785.571419
191 190 0.5286134880 830.512238
192 191 1.1708664900 876.656821
193 192 1.8807594300 924.018452
194 193 2.4199951400 972.610416
195 194 2.5111568000 1022.445998
196 195 4.7797643200 1073.538483
197 196 6.1432360700 1125.901156
198 197 6.4256431400 1179.547301
199 198 4.8200699600 1234.490203
200 199 4.6338606700 1290.743147
and the solution verification shows
> with(res, abs(sum(Val) - sum(Val_est)) < 1e-3 & Val_est[100] < 2)
[1] TRUE