Page 1 of 1
Parameter Estimation by MLE
Posted: Mon Feb 18, 2013 2:23 am
by rbelhach95
Dear all,
Please see the attached code about mle parameter estimation. My model has eight parameters. The problem I am having is with global convergence. it doesn't seem to converge to an absolute maximum, which is attained by my distribution function, even after changing the initial values for the Maximize instruction. I have tested my results using Monte Carlo simulation and computed the bias and mse and they turn out to be very large. Any feedback is very much appreciated.
Rashid
Re: Parameter Estimation by MLE
Posted: Mon Feb 18, 2013 10:50 am
by TomDoan
I'm a bit confused about what you're trying to do. You have an outer
K loop, which looks like that's looping over experiments, but what's the
I loop inside that supposed to do? As written it doesn't do anything but define and redefine a set of FRMLs. Since the MAXIMIZE is after the end of that loop, only the final value of I will matter.
Both for readability and also to make it easier to spot errors, you should pull anything that doesn't change with a trip through the loop outside of it. For instance, this whole block
Code: Select all
compute lambda01s = (muTs-muxs)/sqrt(sigmaT1s^2-((sigmax1s)^(-2)*sigmaTx1s^2))
compute lambda02s = (muTs-muxs)/sqrt(sigmaT2s^2-((sigmax2s)^(-2)*sigmaTx2s^2))
compute lambda11s = -sigmax1s*(1-(sigmax1s)^(-2)*sigmaTx1s)/sqrt(sigmaT1s^2-((sigmax1s)^(-2)*(sigmaTx1s)^2))
compute lambda12s = -sigmax2s*(1-(sigmax2s)^(-2)*sigmaTx2s)/sqrt(sigmaT2s^2-((sigmax2s)^(-2)*(sigmaTx2s)^2))
compute as = -lambda01s/sqrt(1+lambda11s^2)
compute bs = -lambda02s/sqrt(1+lambda12s^2)
compute cs = 1/sqrt(1+lambda11s^2)
compute ds = 1/sqrt(1+lambda12s^2)
compute g3s = (muTs-muxs)/sqrt(sigmaT1s^2+sigmax1s^2+(-2)*sigmaTx1s)
compute g4s = (muxs-muTs)/sqrt(sigmaT2s^2+sigmax2s^2+(-2)*sigmaTx2s)
compute ps = %cdf(g3s)/(%cdf(g3s)+%cdf(g4s))
appears to be fixed across experiments, so should be pulled outside the K loop.
Re: Parameter Estimation by MLE
Posted: Mon Feb 18, 2013 12:52 pm
by rbelhach95
Dear TomDoan,
Thank you for your reply and feedback. Sorry, I am still learning RATS. The K loop is for the Monte Carlo simulation. This is the way I understand my code: the first sets of EWISE (before the I loop) will generate the actual data (50 observations) then I loop will use these 50 obs. to compute the maximum using the MAXIMIZE instruction, now I understand that I should end the I loop after MAXIMIZE, since as you rightly put it MAXIMIZE would read only the last input, but that would generate a maximum for each i in the loop, am I correct? I would prefer using Ewise in all my FRMLS but I don't know how to do it, is the character & inside X(&i) important in this case, I have included it inside the FRMLS instruction but I don't know if it is necessary. As you also may have noticed, I stored all my estimates, bias and mse using SET .....1 N = 0, I have tried to stored them using matrix notation %zeros(A,n) and selecting the kth row using %xrow(A,n) but it doesn't seem to work.
Thank you
Rachid
Re: Parameter Estimation by MLE
Posted: Tue Feb 19, 2013 7:29 pm
by TomDoan
Somewhat ambititious for a newbie.
MAXIMIZE has an implied loop over the data so the FRML just needs to evaluate a single data point (at T, not I). What you want is the following:
Code: Select all
compute M = 50
compute N = 1000
all M
dec vect v(M) simul1(M) simul2(M) u1(M) u2(M) v1(M) v2(M) j1(M) j2(M) j3(N) s(M) x(M)
dec real mux muT sigmax1 sigmaT1 sigmaTx1 sigmax2 sigmaT2 sigmaTx2
set mux_hat 1 N = 0; set muT_hat 1 N = 0; set sigmax1_hat 1 N = 0; set sigmaT1_hat 1 N = 0; set sigmaTx1_hat 1 N = 0
set sigmax2_hat 1 N = 0; set sigmaT2_hat 1 N = 0; set sigmaTx2_hat 1 N = 0
set bias_1 1 N = 0; set bias_2 1 N = 0.; set bias_3 1 N = 0; set bias_4 1 N = 0; set bias_5 1 N = 0
set bias_6 1 N = 0; set bias_7 1 N = 0.; set bias_8 1 N = 0
set mse_1 1 N = 0; set mse_2 1 N = 0; set mse_3 1 N = 0; set mse_4 1 N = 0; set mse_5 1 N = 0
set mse_6 1 N = 0; set mse_7 1 N = 0; set mse_8 1 N = 0
compute muxs = 0.0, muTs = 1.0, sigmax1s = 1.0, sigmaT1s = 0.5, sigmaTx1s = 0.45, sigmax2s = 1.0, sigmaT2s = 1.0, sigmaTx2s = 0.5
compute lambda01s = (muTs-muxs)/sqrt(sigmaT1s^2-((sigmax1s)^(-2)*sigmaTx1s^2))
compute lambda02s = (muTs-muxs)/sqrt(sigmaT2s^2-((sigmax2s)^(-2)*sigmaTx2s^2))
compute lambda11s = -sigmax1s*(1-(sigmax1s)^(-2)*sigmaTx1s)/sqrt(sigmaT1s^2-((sigmax1s)^(-2)*(sigmaTx1s)^2))
compute lambda12s = -sigmax2s*(1-(sigmax2s)^(-2)*sigmaTx2s)/sqrt(sigmaT2s^2-((sigmax2s)^(-2)*(sigmaTx2s)^2))
compute as = -lambda01s/sqrt(1+lambda11s^2)
compute bs = -lambda02s/sqrt(1+lambda12s^2)
compute cs = 1/sqrt(1+lambda11s^2)
compute ds = 1/sqrt(1+lambda12s^2)
compute g3s = (muTs-muxs)/sqrt(sigmaT1s^2+sigmax1s^2+(-2)*sigmaTx1s)
compute g4s = (muxs-muTs)/sqrt(sigmaT2s^2+sigmax2s^2+(-2)*sigmaTx2s)
compute ps = %cdf(g3s)/(%cdf(g3s)+%cdf(g4s))
nonlin mux muT sigmax1 sigmaT1 sigmaTx1 sigmax2 sigmaT2 sigmaTx2
do k = 1,N
ewise v(i)= %ran(1.0)
ewise simul1(i)= %uniform(0,1)
ewise u1(i) = %invnormal(%cdf(as)+((1-%cdf(as))*simul1(i)))
ewise u2(i) = %invnormal(%cdf(bs)*simul1(i))
ewise simul2(i)= %uniform(0,1)
ewise j1(i) = %sign(simul2(i)-ps)
ewise j2(i) = j1(i)+1
ewise j3(i) = %sign(j2(i))
ewise s(i) = cs*(lambda11s*u1(i)+v(i))*(1-j3(i)) + ds*(lambda12s*u2(i)+v(i))*j3(i)
ewise x(i) = (sigmax1s*s(i)+muxs)*(1-j3(i)) + (sigmax2s*s(i)+muxs)*j3(i)
frml y = $
g1 =(1/sigmax1)*%density((x(t)-mux)/sigmax1),$
g21 = -((1-((sigmax1)^(-2)*sigmaTx1))*(x(t)-mux))/sqrt(sigmaT1^2-(sigmax1^(-2)*sigmaTx1^2)),$
g22 = (muT-mux)/sqrt(sigmaT1^2-(sigmax1^(-2)*(sigmaTx1)^2)),$
g2 = %cdf(g21+g22),$
g3 = %cdf((muT-mux)/sqrt(sigmaT1^2+(sigmax1)^2+(-2)*sigmaTx1)),$
g4 = %cdf((mux-muT)/sqrt(sigmaT2^2+(sigmax2)^2 +(-2)*sigmaTx2)),$
h1 = (1/sigmax2)*%density((x(t)-mux)/sigmax2),$
h21 = ((1-((sigmax2)^(-2)*sigmaTx2))*(x(t)-mux))/sqrt((sigmaT2^2-(sigmax2^(-2)*sigmaTx2^2))),$
h22 = (mux-muT)/sqrt(sigmaT2^2-(sigmax2^(-2)*sigmaTx2^2)),$
h2 = %cdf(h21+h22),$
log((g1*g2)+(h1*h2)) - log(g3+g4)
compute mux = 0.0, muT = 1.0, sigmax1 = 1.0, sigmaT1 = 0.5, sigmaTx1 = 0.45, sigmax2 = 1.0, sigmaT2 = 1.0, sigmaTx2 = 0.5
maximize(pmethod=simplex,piters=9,iters = 5000, method=bfgs, cvcrit=0.0001) y
com mux_hat(k) = %beta(1), muT_hat(k) = %beta(2)
com sigmax1_hat(k) = %beta(3), sigmaT1_hat(k) = %beta(4)
com sigmax2_hat(k) = %beta(5), sigmaT2_hat(k) = %beta(6)
com sigmaTx1_hat(k) = %beta(7), sigmaTx2_hat(k) = %beta(8)
com bias_1(k) = muxs - mux_hat(k), bias_2(k) = muTs - muT_hat(k)
com bias_3(k) = sigmax1s - sigmax1_hat(k), bias_4(k) = sigmaT1s - sigmaT1_hat(k)
com bias_5(k) = sigmax2s - sigmax2_hat(k), bias_6(k) = sigmaT2s - sigmaT2_hat(k)
com bias_7(k) = sigmaTx1s - sigmaTx1_hat(k), bias_8(k) = sigmaTx2s - sigmaTx2_hat(k)
com mse_1(k) = (bias_1(k))^2, mse_2(k) = (bias_2(k))^2
com mse_3(k) = (bias_3(k))^2, mse_4(k) = (bias_4(k))^2, mse_5(k) = (bias_5(k))^2
com mse_6(k) = (bias_6(k))^2, mse_7(k) = (bias_7(k))^2, mse_8(k) = (bias_8(k))^2
end do k
table / bias_1 bias_2 bias_3 bias_4 bias_5 bias_6 bias_7 bias_8
table / mse_1 mse_2 mse_3 mse_4 mse_5 mse_6 mse_7 mse_8
You might want to double-check that the formula is what you want. The optimization has rather odd behavior, as BFGS has a great deal of difficulty improving on even a small number of simplex iterations; almost like it's bumping into a discontinuity.
Re: Parameter Estimation by MLE
Posted: Wed Feb 20, 2013 4:34 pm
by rbelhach95
Thanks for your reply, very much appreciated. It is still not working as it should. My distribution is not just continuous but infinitely differentiable too. my mse and bias are rather large. Could you check to see the way I stored the eight estimates, is it right? I would prefer to use matrices but I have tried it once like I explained in my last email but to no avail. Thank you
Rachid
Re: Parameter Estimation by MLE
Posted: Thu Feb 21, 2013 4:03 am
by rbelhach95
Dear TomDoan,
I have noticed that my post MLE Parameters Estimation has disappeared from the Examples and Sample code, please advise!
Rachid