Beefy Boxes and Bandwidth Generously Provided by pair Networks
go ahead... be a heretic
 
PerlMonks  

comment on

( [id://3333]=superdoc: print w/replies, xml ) Need Help??

grim_master:

There's too much domain-specific jargon for me to know what the heck the program is supposed to compute. However, I went through your program and removed unnecessary operations, simplifying as much as I could. (I tend to do this to figure out what programs are doing.)

In the process of playing with your program, I noticed a few things:

  • You only go through your loop twice. Glancing at the paper, I'm expecting you'll need a larger set of iterations for it to do anything useful.
  • Structurally, your program is dangerous. You use global variables which you modify in and out of your subroutine. It's very easy to get confused about what's happening, and you could be having trouble with "action at a distance". (Meaning, you may be altering a value when you're not expecting to do so.)
  • Something that will help you keep things clear is to declare variables:
    • In the smallest scope you can, and
    • only when and where you need to.
    For example: You're declaring a global $result variable, whose sole purpose is to sit inside your lambda subroutine. In there, it simply receives the sum of two values, and the next statement returns the value. I simply removed the $result variable and returned the result directly. (Later, I extracted $lambda0 from the subroutine to remove the "action at a distance".
  • Your lack of indentation makes the program impossible to read, so the first thing I had to do was fix the indentation.
  • You wrote your equations in a rather obtuse style, making them harder to look at. I trimmed them a bit.
  • In your for loop, you have an if statement to separate the initial conditions setup from the actual code. It's cleaner to simply set up the initial conditions before the loop, and adjust your starting index accordingly. That way, you don't have to look "around" the if statement each time through the loop.

I could continue, but I don't really have the time. (Nor the inclination, as I don't really understand the paper.) So the code is still a bit of a mess, but here's what I did with it. I definitely oversimplified it (by removing the $i variable, for example). You'll probably want to put that back in when you're fixing it. Now that I've properly disclaimed it, here's what I simplified it to:

#!/usr/bin/perl use strict; use warnings; use diagnostics; #Force the same set of random numbers for each run until it works... #use Math::Random::MT qw(rand); srand 0; open (F1,'>>event_times_new2.mcm'); open (F2,'>>lambda_star_new2.mcm'); my ($lambda0, $alpha, $beta, $lambda_star); my $time=8000000; my $counterfails=0; my $num=1; $lambda0=1.2; $alpha=0.8; $beta=0.9; #ALGO it self my (@d, @s, @u, @event_times); $lambda_star = $lambda0; $d[0] = 0; $u[0] = rand($num); $s[0] = -log($u[0]) / $lambda_star; die if $s[0] >= $time; $event_times[0] = $s[0]; print F1 "$event_times[0]\n"; print F2 "$lambda_star\n"; my $valuelambda = lambda($event_times[0]) + $lambda0; $lambda_star=$valuelambda+$alpha; START: $u[1]=rand($num); $s[1] = $s[0] - log($u[1])/$lambda_star; die if $s[1] > $time; $d[1]=rand($num) ; my $rejectionlambda=lambda($s[1]); if ($d[1] < $rejectionlambda/$lambda_star) { $event_times[1]=$s[1]; print F1 "$event_times[1]\n"; print F2 "$lambda_star\n"; }else{ $counterfails+=1; $lambda_star=$rejectionlambda; goto START; } print F1 "The number of event-times that failed the A/R procedure is $ +counterfails\n"; close F1; close F2; ################################################ #subroutine that computes intensity function sub lambda { my $t = $_[0]; my $term; my $sum=0; for (my $i=0; $i<@event_times; $i++) { $term = $alpha*exp( -$beta * ($t - $event_times[$i]) ); $sum += $term; } $sum -= $alpha if $term == $alpha; return $sum; }

NOTE: One trick I used when modifying your program is to use the same sequence of random numbers for each run. Since I'm just learning about your program, I didn't want to change the values returned. So by using the same sequence of random numbers each time, I could look at the results after each change and see if they changed. If so, then I made a mistake, so I back out the change.

If you're curious about any particular change(s) I made, just ask. I'll be happy to tell you why I did each one.

...roboticus

When your only tool is a hammer, all problems look like your thumb.


In reply to Re: Simulating uni-dimensional Hawkes processes by roboticus
in thread Simulating uni-dimensional Hawkes processes by glrm_master

Title:
Use:  <p> text here (a paragraph) </p>
and:  <code> code here </code>
to format your post; it's "PerlMonks-approved HTML":



  • Are you posting in the right place? Check out Where do I post X? to know for sure.
  • Posts may use any of the Perl Monks Approved HTML tags. Currently these include the following:
    <code> <a> <b> <big> <blockquote> <br /> <dd> <dl> <dt> <em> <font> <h1> <h2> <h3> <h4> <h5> <h6> <hr /> <i> <li> <nbsp> <ol> <p> <small> <strike> <strong> <sub> <sup> <table> <td> <th> <tr> <tt> <u> <ul>
  • Snippets of code should be wrapped in <code> tags not <pre> tags. In fact, <pre> tags should generally be avoided. If they must be used, extreme care should be taken to ensure that their contents do not have long lines (<70 chars), in order to prevent horizontal scrolling (and possible janitor intervention).
  • Want more info? How to link or How to display code and escape characters are good places to start.
Log In?
Username:
Password:

What's my password?
Create A New User
Domain Nodelet?
Chatterbox?
and the web crawler heard nothing...

How do I use this?Last hourOther CB clients
Other Users?
Others making s'mores by the fire in the courtyard of the Monastery: (3)
As of 2024-04-26 02:16 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    No recent polls found