function part2sim2pathssin(K,L,N,sigma,numobs)
  
  %P is here assumed known 
  
  angl = pi/4;
  d = 1;
  lambda = 10;
  
  I2 = lambda*log( (1+sin(angl))/(1-sin(angl)) ) / ( 4*d*(angl^2) );
  I3 = (lambda^2)*tan(angl)/(4*(d^2)*(angl^3));
  
  P1 = (1/3)*(0.5+1+1.5);
  P2 = (1/3)*(0.5^2 + 1 + 1.5^2);
  P3 = (1/3)*(0.5^3 + 1 + 1.5^3);
  
  W1 = [];
  W2 = [];
  W3 = [];
  for (Ltry=1:max(K,N))
    c1 = N/K;
    c2 = Ltry/N;
    c3 = Ltry/K;
  
    T1   = P1;
    T2   = P2 + c2*I2*P1^2;
    T3   = P3 + 3*c2*I2*P1*P2 + c2^2*I3*P1^3;
    T11  = P1^2;
    T111 = P1^3;
    T12  = P1*P2 + c2*I2*P1^3;
  
    S1 = T1;
    S2 = T2 + c3*T11;
    S3 = (1 + 1/(K^2))*T3 + 3*c3*T12 + c3^2*T111;
    S11 = T11 + (1/(K*Ltry))*T2;
  
    gamma1 = c2*S1;
    gamma2 = c2*S2;
    gamma3 = c2*S3;
    gamma11 = c2^2*S11;
  
    W1(Ltry) = gamma1 + sigma^2;
    W2(Ltry) = gamma2 + 2*sigma^2*(1+c1)*gamma1 + sigma^4*(1 + c1);
    W3(Ltry) = gamma3 + 3*sigma^2*(1+c1)*gamma2 + 3*sigma^2*c1*gamma11 + 3*sigma^4*( c1^2+3*c1+1+1/(K^2) )*gamma1 + sigma^6*( c1^2+3*c1+1+1/(K^2) );
  end
  
  P = zeros(L,L);
  step = floor(L/3);
  for (k=1:step)
    P(k,k) = 0.5;
    P(k+step,k+step) = 1;
    P(k+2*step,k+2*step) = 1.5;
  end
  
  closestL = [];
  for (obsnr=1:numobs)
    obsnr
    samplecovsum1 = 0;
    samplecovsum2 = 0;
    samplecovsum3 = 0;
    for (thisobs=1:obsnr)
      V = zeros(N,L);
      for (k=1:L) 
        theta = 2*angl*(rand(1)-0.5);
        V(:,k) = ((exp(2*pi*j*sin(theta)*d/lambda)).^(0:(N-1)))';
      end
      V = (1/sqrt(N))*V;
      alpha = randn(L,K);
      nmatrix = sigma*randn(N,K);
      Y = V*sqrt(P)*alpha + nmatrix;
      W = (1/K)*Y*Y';%sample covariance matrix sample
      samplecovsum1 = samplecovsum1 + (1/N)*trace(W);
      samplecovsum2 = samplecovsum2 + (1/N)*trace(W^2);
      samplecovsum3 = samplecovsum3 + (1/N)*trace(W^3);
    end
    
    closestL(obsnr) = 1;
    leastsse = ((1/obsnr)*samplecovsum1 - W1(1))^2 + ((1/obsnr)*samplecovsum2 - W2(1))^2 + ((1/obsnr)*samplecovsum3 - W3(1))^2;
    for (tryL=2:max(K,N))
      sse = ((1/obsnr)*samplecovsum1 - W1(tryL))^2 + ((1/obsnr)*samplecovsum2 - W2(tryL))^2 + ((1/obsnr)*samplecovsum3 - W3(tryL))^2;
      if (sse < leastsse)
        leastsse = sse;
        closestL(obsnr) = tryL;
      end
    end
  end
  
  plot(1:numobs,closestL,'ko',1:numobs,L*ones(1,numobs),'k-');
  xlabel('Number of observations')
  ylabel('L');
  axis([0 numobs 0 2*L]);
  legend('Estimate of L','Actual value of L');