Help with basic code

I attached the pdf its the last page 9 I meant to code that in MATLAB, I believe even the code is basic but I am not sure why is not working. I tried to change the code several time but no luck.

http://www.stat.sdu.dk/matstat/yuri-st505/w5slides1.pdf

from the Paper:

(1) t = 0, I = 0.
(2) Generate a random number U ~ U(0, 1).
(3) t = t - (1/Lam)*ln(U). If t > T then stop.
(4) Generate a random number U ~ U(0, 1).
(5) If U <= Lam(t)/Lam, set I = I + 1, S(I) = t.
(6) Go to step 2.

Output:
I: the number of events at time T ,
S(1), . . . , S(I): the event times.
2. Direct generation of successive event times

function nonhomogeneous (i,s)
i = 0;
t = 010;
T = length(cos(t));
s = [0,t];
u = rand (1,1);
t = t - (1 /cos(t)*log(u));
while t <= T
    u2 = rand (1,1);
    if ( u2 <= cos(t) / max(cos(t)))
        hold on
        plot(t,cos(t))
        i = i+1;        
        s(i)=t;
        u = rand (1,1);
        t = t - (1 /cos(t)*log(u));

      end
end

Answer 1

[Andrew Newell]

Here is some code that solves the problem and also illustrates an important programming principle: when you create a function, input the variables that you want to experiment with.

Here is a function that outputs a nonhomogeneous random sequence of times.

function S = nonhomogeneous(lambda0,lambda,T)
t = 0;
I = 0;
S = [];

 u = rand;
t = t - log(u)/lambda0;
 while t <= T
    u = rand;
    if (u <= lambda(t)/lambda0)
        I = I+1;
        S(I) = t;
    end
    u = rand;
    t = t - log(u)/lambda0;
end

Save this in a file nonhomogenous.m. Here is a script to run it:

lambda0 = 50; % Needs to be much larger than T to get several numbers
T = 1;

lambda = @(x) lambda0*cos(x); % This is lambda0(t)/lambda0
S = nonhomogeneous(lambda0,lambda,T);
subplot 121
plot(S,lambda(S),'.')
xlabel('t')
ylabel('lambda(t)')

lambda = @(x) lambda0*sin(x);
S = nonhomogeneous(lambda0,lambda,T);
subplot 122
plot(S,lambda(S),'.')
xlabel('t')

I'm not sure what you mean by "spikes". Maybe you really want to plot something like

hist(S)

[Susan]

Thanks very much I really appreciate it.. the only question I have is theat after the if statement u is <= f(t)/ then we should divide it by the lambda ( and I believe the lambda is the maximum value of lambda) but in your code you didn't do that, other than this point thanks for your help I will debug it to understand it exactly

if (u <= f(t))

I = I+1;

S(I) = t;

end

[Andrew Newell]

f(t) is lambda(t)/lambda. I should have written it this way for clarity, so I have edited the above code. Note that lambda0 is the maximum value for lambda in the code.

[Susan]

Thanks very much :)

Answer 2

[Andrew Newell]

@Susan, thank you for formatting the code nicely. I wish more people would in their questions!

I notice your function doesn't have any outputs and your inputs are redundant. If you change the top line to

function [s,i] = nonhomogeneous

do you get what you want?

EDIT: I see several problems:

t = 010; % t=10

instead of

t = 0;

Also,

T = length(cos(t));

gives T=1 (is that what you're after?)

s = [0,t];

initializes your array s to [0 10] instead of [].

Yet another: there is no lambda. Instead, you're using

max(cos(t))

which is only the maximum for the current values of t, not 1 as it should be.

[Susan]

Thanks very much Ive changed it accordingly but still I dont get what I am expecting.. I am meant to make the lambda function e.g sin or cos and then to have a graph at the end with that function and the higher the curve of the function the more spikes and vice versa..

the plot is not working at all.

Answer 3

[Jan Jensen]

I've read the presentation, and it was quite interesting. Not sure I understand all the implications, though.

Please post your definition of the lambda function as well as the code.

A good starting point could also be to implement the simulation of Poisson process with constant decay (lambda = constant). The one on slide 5. The dimensions of the output variable S should be the same. If the plot of S for lambda=constant makes sense you have made a great step forward in debugging your code, including your plotting commands :)

[Susan]

Hey Jan..

I tried it with constant and it is workingm I just changed few things but its working finally :)