Skip to content

Commit

Permalink
Deploy to GitHub Pages: a3d30df
Browse files Browse the repository at this point in the history
  • Loading branch information
@kaskr
kaskr committed Nov 4, 2023
1 parent a2fff21 commit 4ac585c
Showing 1 changed file with 1 addition and 1 deletion.
2 changes: 1 addition & 1 deletion tweedie_8cpp_source.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,6 @@
<div class="title">tweedie.cpp</div> </div>
</div><!--header-->
<div class="contents">
<div class="fragment"><div class="line"><a name="l00001"></a><span class="lineno"> 1</span>&#160;<span class="comment">// Re-structured version of tweedie density function from &#39;cplm&#39; package.</span></div><div class="line"><a name="l00002"></a><span class="lineno"> 2</span>&#160;</div><div class="line"><a name="l00003"></a><span class="lineno"> 3</span>&#160;<span class="comment">/* the threshold used in finding the bounds of the series */</span></div><div class="line"><a name="l00004"></a><span class="lineno"> 4</span>&#160;<span class="preprocessor">#define TWEEDIE_DROP 37.0</span></div><div class="line"><a name="l00005"></a><span class="lineno"> 5</span>&#160;<span class="comment">/* the loop increment used in finding the bounds of the series */</span></div><div class="line"><a name="l00006"></a><span class="lineno"> 6</span>&#160;<span class="preprocessor">#define TWEEDIE_INCRE 5</span></div><div class="line"><a name="l00007"></a><span class="lineno"> 7</span>&#160;<span class="comment">/* the max number of terms allowed in the finite sum approximation*/</span></div><div class="line"><a name="l00008"></a><span class="lineno"> 8</span>&#160;<span class="preprocessor">#define TWEEDIE_NTERM 20000</span></div><div class="line"><a name="l00009"></a><span class="lineno"> 9</span>&#160;</div><div class="line"><a name="l00018"></a><span class="lineno"> 18</span>&#160;<span class="keyword">template</span>&lt;<span class="keyword">class</span> Float&gt;</div><div class="line"><a name="l00019"></a><span class="lineno"> 19</span>&#160;Float tweedie_logW(Float y, Float phi, Float p){</div><div class="line"><a name="l00020"></a><span class="lineno"> 20</span>&#160; <span class="keywordtype">bool</span> ok = (0 &lt; y) &amp;&amp; (0 &lt; phi) &amp;&amp; (1 &lt; p) &amp;&amp; (p &lt; 2);</div><div class="line"><a name="l00021"></a><span class="lineno"> 21</span>&#160; <span class="keywordflow">if</span> (!ok) <span class="keywordflow">return</span> NAN;</div><div class="line"><a name="l00022"></a><span class="lineno"> 22</span>&#160;</div><div class="line"><a name="l00023"></a><span class="lineno"> 23</span>&#160; Float p1 = p - 1.0, p2 = 2.0 - p;</div><div class="line"><a name="l00024"></a><span class="lineno"> 24</span>&#160; Float a = - p2 / p1, a1 = 1.0 / p1;</div><div class="line"><a name="l00025"></a><span class="lineno"> 25</span>&#160; Float cc, w, sum_ww = 0.0;</div><div class="line"><a name="l00026"></a><span class="lineno"> 26</span>&#160; <span class="keywordtype">double</span> ww_max = -INFINITY ;</div><div class="line"><a name="l00027"></a><span class="lineno"> 27</span>&#160; <span class="keywordtype">double</span> j;</div><div class="line"><a name="l00028"></a><span class="lineno"> 28</span>&#160;</div><div class="line"><a name="l00029"></a><span class="lineno"> 29</span>&#160; <span class="comment">/* only need the lower bound and the # terms to be stored */</span></div><div class="line"><a name="l00030"></a><span class="lineno"> 30</span>&#160; <span class="keywordtype">int</span> jh, jl, jd;</div><div class="line"><a name="l00031"></a><span class="lineno"> 31</span>&#160; <span class="keywordtype">double</span> jmax = 0;</div><div class="line"><a name="l00032"></a><span class="lineno"> 32</span>&#160; Float logz = 0;</div><div class="line"><a name="l00033"></a><span class="lineno"> 33</span>&#160;</div><div class="line"><a name="l00034"></a><span class="lineno"> 34</span>&#160; <span class="comment">/* compute jmax for the given y &gt; 0*/</span></div><div class="line"><a name="l00035"></a><span class="lineno"> 35</span>&#160; cc = a * log(p1) - log(p2);</div><div class="line"><a name="l00036"></a><span class="lineno"> 36</span>&#160; jmax = asDouble( fmax2(1.0, pow(y, p2) / (phi * p2)) );</div><div class="line"><a name="l00037"></a><span class="lineno"> 37</span>&#160; logz = - a * log(y) - a1 * log(phi) + cc;</div><div class="line"><a name="l00038"></a><span class="lineno"> 38</span>&#160;</div><div class="line"><a name="l00039"></a><span class="lineno"> 39</span>&#160; <span class="comment">/* find bounds in the summation */</span></div><div class="line"><a name="l00040"></a><span class="lineno"> 40</span>&#160; <span class="comment">/* locate upper bound */</span></div><div class="line"><a name="l00041"></a><span class="lineno"> 41</span>&#160; cc = logz + a1 + a * log(-a);</div><div class="line"><a name="l00042"></a><span class="lineno"> 42</span>&#160; j = jmax ;</div><div class="line"><a name="l00043"></a><span class="lineno"> 43</span>&#160; w = a1 * j ;</div><div class="line"><a name="l00044"></a><span class="lineno"> 44</span>&#160; <span class="keywordflow">while</span> (1) {</div><div class="line"><a name="l00045"></a><span class="lineno"> 45</span>&#160; j += TWEEDIE_INCRE ;</div><div class="line"><a name="l00046"></a><span class="lineno"> 46</span>&#160; <span class="keywordflow">if</span> (j * (cc - a1 * log(j)) &lt; (w - TWEEDIE_DROP))</div><div class="line"><a name="l00047"></a><span class="lineno"> 47</span>&#160; break ;</div><div class="line"><a name="l00048"></a><span class="lineno"> 48</span>&#160; }</div><div class="line"><a name="l00049"></a><span class="lineno"> 49</span>&#160; jh = ceil(j);</div><div class="line"><a name="l00050"></a><span class="lineno"> 50</span>&#160; <span class="comment">/* locate lower bound */</span></div><div class="line"><a name="l00051"></a><span class="lineno"> 51</span>&#160; j = jmax;</div><div class="line"><a name="l00052"></a><span class="lineno"> 52</span>&#160; <span class="keywordflow">while</span> (1) {</div><div class="line"><a name="l00053"></a><span class="lineno"> 53</span>&#160; j -= TWEEDIE_INCRE ;</div><div class="line"><a name="l00054"></a><span class="lineno"> 54</span>&#160; <span class="keywordflow">if</span> (j &lt; 1 || j * (cc - a1 * log(j)) &lt; w - TWEEDIE_DROP)</div><div class="line"><a name="l00055"></a><span class="lineno"> 55</span>&#160; break ;</div><div class="line"><a name="l00056"></a><span class="lineno"> 56</span>&#160; }</div><div class="line"><a name="l00057"></a><span class="lineno"> 57</span>&#160; jl = imax2(1, floor(j)) ;</div><div class="line"><a name="l00058"></a><span class="lineno"> 58</span>&#160; jd = jh - jl + 1;</div><div class="line"><a name="l00059"></a><span class="lineno"> 59</span>&#160;</div><div class="line"><a name="l00060"></a><span class="lineno"> 60</span>&#160; <span class="comment">/* set limit for # terms in the sum */</span></div><div class="line"><a name="l00061"></a><span class="lineno"> 61</span>&#160; <span class="keywordtype">int</span> nterms = imin2(jd, TWEEDIE_NTERM) ;</div><div class="line"><a name="l00062"></a><span class="lineno"> 62</span>&#160; std::vector&lt;Float&gt; ww(nterms);</div><div class="line"><a name="l00063"></a><span class="lineno"> 63</span>&#160; <span class="comment">/* evaluate series using the finite sum*/</span></div><div class="line"><a name="l00064"></a><span class="lineno"> 64</span>&#160; <span class="comment">/* y &gt; 0 */</span></div><div class="line"><a name="l00065"></a><span class="lineno"> 65</span>&#160; sum_ww = 0.0 ;</div><div class="line"><a name="l00066"></a><span class="lineno"> 66</span>&#160; <span class="keywordtype">int</span> iterm = imin2(jd, nterms) ;</div><div class="line"><a name="l00067"></a><span class="lineno"> 67</span>&#160; <span class="keywordflow">for</span> (<span class="keywordtype">int</span> k = 0; k &lt; iterm; k++) {</div><div class="line"><a name="l00068"></a><span class="lineno"> 68</span>&#160; j = k + jl ;</div><div class="line"><a name="l00069"></a><span class="lineno"> 69</span>&#160; ww[k] = j * logz - <a class="code" href="group__special__functions.html#ga177b023eb871927fcfc90259e9f00382">lgamma</a>(1 + j) - <a class="code" href="group__special__functions.html#ga177b023eb871927fcfc90259e9f00382">lgamma</a>(-a * j);</div><div class="line"><a name="l00070"></a><span class="lineno"> 70</span>&#160; ww_max = fmax2( ww_max, asDouble(ww[k]) );</div><div class="line"><a name="l00071"></a><span class="lineno"> 71</span>&#160; }</div><div class="line"><a name="l00072"></a><span class="lineno"> 72</span>&#160; <span class="keywordflow">for</span> (<span class="keywordtype">int</span> k = 0; k &lt; iterm; k++)</div><div class="line"><a name="l00073"></a><span class="lineno"> 73</span>&#160; sum_ww += exp(ww[k] - ww_max);</div><div class="line"><a name="l00074"></a><span class="lineno"> 74</span>&#160; Float ans = log(sum_ww) + ww_max ;</div><div class="line"><a name="l00075"></a><span class="lineno"> 75</span>&#160;</div><div class="line"><a name="l00076"></a><span class="lineno"> 76</span>&#160; <span class="keywordflow">return</span> ans;</div><div class="line"><a name="l00077"></a><span class="lineno"> 77</span>&#160;}</div><div class="ttc" id="group__special__functions_html_ga177b023eb871927fcfc90259e9f00382"><div class="ttname"><a href="group__special__functions.html#ga177b023eb871927fcfc90259e9f00382">lgamma</a></div><div class="ttdeci">Type lgamma(Type x)</div><div class="ttdoc">Logarithm of gamma function (following R argument convention). </div><div class="ttdef"><b>Definition:</b> <a href="lgamma_8hpp_source.html#l00012">lgamma.hpp:12</a></div></div>
<div class="fragment"><div class="line"><a name="l00001"></a><span class="lineno"> 1</span>&#160;<span class="comment">// Re-structured version of tweedie density function from &#39;cplm&#39; package.</span></div><div class="line"><a name="l00002"></a><span class="lineno"> 2</span>&#160;</div><div class="line"><a name="l00003"></a><span class="lineno"> 3</span>&#160;<span class="comment">/* the threshold used in finding the bounds of the series */</span></div><div class="line"><a name="l00004"></a><span class="lineno"> 4</span>&#160;<span class="preprocessor">#define TWEEDIE_DROP 37.0</span></div><div class="line"><a name="l00005"></a><span class="lineno"> 5</span>&#160;<span class="comment">/* the loop increment used in finding the bounds of the series */</span></div><div class="line"><a name="l00006"></a><span class="lineno"> 6</span>&#160;<span class="preprocessor">#define TWEEDIE_INCRE 5</span></div><div class="line"><a name="l00007"></a><span class="lineno"> 7</span>&#160;<span class="comment">/* the max number of terms allowed in the finite sum approximation*/</span></div><div class="line"><a name="l00008"></a><span class="lineno"> 8</span>&#160;<span class="preprocessor">#define TWEEDIE_NTERM 20000</span></div><div class="line"><a name="l00009"></a><span class="lineno"> 9</span>&#160;</div><div class="line"><a name="l00018"></a><span class="lineno"> 18</span>&#160;<span class="keyword">template</span>&lt;<span class="keyword">class</span> Float&gt;</div><div class="line"><a name="l00019"></a><span class="lineno"> 19</span>&#160;Float tweedie_logW(Float y, Float phi, Float p){</div><div class="line"><a name="l00020"></a><span class="lineno"> 20</span>&#160; <span class="keywordtype">bool</span> ok = (0 &lt; y) &amp;&amp; (0 &lt; phi) &amp;&amp; (1 &lt; p) &amp;&amp; (p &lt; 2);</div><div class="line"><a name="l00021"></a><span class="lineno"> 21</span>&#160; <span class="keywordflow">if</span> (!ok) <span class="keywordflow">return</span> NAN;</div><div class="line"><a name="l00022"></a><span class="lineno"> 22</span>&#160;</div><div class="line"><a name="l00023"></a><span class="lineno"> 23</span>&#160; Float p1 = p - 1.0, p2 = 2.0 - p;</div><div class="line"><a name="l00024"></a><span class="lineno"> 24</span>&#160; Float a = - p2 / p1, a1 = 1.0 / p1;</div><div class="line"><a name="l00025"></a><span class="lineno"> 25</span>&#160; Float cc, w, sum_ww = 0.0;</div><div class="line"><a name="l00026"></a><span class="lineno"> 26</span>&#160; <span class="keywordtype">double</span> ww_max = -INFINITY ;</div><div class="line"><a name="l00027"></a><span class="lineno"> 27</span>&#160; <span class="keywordtype">double</span> j;</div><div class="line"><a name="l00028"></a><span class="lineno"> 28</span>&#160;</div><div class="line"><a name="l00029"></a><span class="lineno"> 29</span>&#160; <span class="comment">/* only need the lower bound and the # terms to be stored */</span></div><div class="line"><a name="l00030"></a><span class="lineno"> 30</span>&#160; <span class="keywordtype">double</span> jh, jl, jd;</div><div class="line"><a name="l00031"></a><span class="lineno"> 31</span>&#160; <span class="keywordtype">double</span> jmax = 0;</div><div class="line"><a name="l00032"></a><span class="lineno"> 32</span>&#160; Float logz = 0;</div><div class="line"><a name="l00033"></a><span class="lineno"> 33</span>&#160;</div><div class="line"><a name="l00034"></a><span class="lineno"> 34</span>&#160; <span class="comment">/* compute jmax for the given y &gt; 0*/</span></div><div class="line"><a name="l00035"></a><span class="lineno"> 35</span>&#160; cc = a * log(p1) - log(p2);</div><div class="line"><a name="l00036"></a><span class="lineno"> 36</span>&#160; jmax = asDouble( fmax2(1.0, pow(y, p2) / (phi * p2)) );</div><div class="line"><a name="l00037"></a><span class="lineno"> 37</span>&#160; logz = - a * log(y) - a1 * log(phi) + cc;</div><div class="line"><a name="l00038"></a><span class="lineno"> 38</span>&#160;</div><div class="line"><a name="l00039"></a><span class="lineno"> 39</span>&#160; <span class="comment">/* find bounds in the summation */</span></div><div class="line"><a name="l00040"></a><span class="lineno"> 40</span>&#160; <span class="comment">/* locate upper bound */</span></div><div class="line"><a name="l00041"></a><span class="lineno"> 41</span>&#160; cc = logz + a1 + a * log(-a);</div><div class="line"><a name="l00042"></a><span class="lineno"> 42</span>&#160; j = jmax ;</div><div class="line"><a name="l00043"></a><span class="lineno"> 43</span>&#160; w = a1 * j ;</div><div class="line"><a name="l00044"></a><span class="lineno"> 44</span>&#160; <span class="keywordflow">while</span> (1) {</div><div class="line"><a name="l00045"></a><span class="lineno"> 45</span>&#160; j += TWEEDIE_INCRE ;</div><div class="line"><a name="l00046"></a><span class="lineno"> 46</span>&#160; <span class="keywordflow">if</span> (j * (cc - a1 * log(j)) &lt; (w - TWEEDIE_DROP))</div><div class="line"><a name="l00047"></a><span class="lineno"> 47</span>&#160; break ;</div><div class="line"><a name="l00048"></a><span class="lineno"> 48</span>&#160; }</div><div class="line"><a name="l00049"></a><span class="lineno"> 49</span>&#160; jh = ceil(j);</div><div class="line"><a name="l00050"></a><span class="lineno"> 50</span>&#160; <span class="comment">/* locate lower bound */</span></div><div class="line"><a name="l00051"></a><span class="lineno"> 51</span>&#160; j = jmax;</div><div class="line"><a name="l00052"></a><span class="lineno"> 52</span>&#160; <span class="keywordflow">while</span> (1) {</div><div class="line"><a name="l00053"></a><span class="lineno"> 53</span>&#160; j -= TWEEDIE_INCRE ;</div><div class="line"><a name="l00054"></a><span class="lineno"> 54</span>&#160; <span class="keywordflow">if</span> (j &lt; 1 || j * (cc - a1 * log(j)) &lt; w - TWEEDIE_DROP)</div><div class="line"><a name="l00055"></a><span class="lineno"> 55</span>&#160; break ;</div><div class="line"><a name="l00056"></a><span class="lineno"> 56</span>&#160; }</div><div class="line"><a name="l00057"></a><span class="lineno"> 57</span>&#160; jl = std::fmax(1., floor(j)) ;</div><div class="line"><a name="l00058"></a><span class="lineno"> 58</span>&#160; jd = jh - jl + 1;</div><div class="line"><a name="l00059"></a><span class="lineno"> 59</span>&#160;</div><div class="line"><a name="l00060"></a><span class="lineno"> 60</span>&#160; <span class="comment">/* set limit for # terms in the sum */</span></div><div class="line"><a name="l00061"></a><span class="lineno"> 61</span>&#160; <span class="keywordtype">int</span> nterms = (int) std::fmin(jd, TWEEDIE_NTERM) ;</div><div class="line"><a name="l00062"></a><span class="lineno"> 62</span>&#160; std::vector&lt;Float&gt; ww(nterms);</div><div class="line"><a name="l00063"></a><span class="lineno"> 63</span>&#160; <span class="comment">/* evaluate series using the finite sum*/</span></div><div class="line"><a name="l00064"></a><span class="lineno"> 64</span>&#160; <span class="comment">/* y &gt; 0 */</span></div><div class="line"><a name="l00065"></a><span class="lineno"> 65</span>&#160; sum_ww = 0.0 ;</div><div class="line"><a name="l00066"></a><span class="lineno"> 66</span>&#160; <span class="keywordtype">int</span> iterm = (int) std::fmin(jd, nterms) ;</div><div class="line"><a name="l00067"></a><span class="lineno"> 67</span>&#160; <span class="keywordflow">for</span> (<span class="keywordtype">int</span> k = 0; k &lt; iterm; k++) {</div><div class="line"><a name="l00068"></a><span class="lineno"> 68</span>&#160; j = k + jl ;</div><div class="line"><a name="l00069"></a><span class="lineno"> 69</span>&#160; ww[k] = j * logz - <a class="code" href="group__special__functions.html#ga177b023eb871927fcfc90259e9f00382">lgamma</a>(1 + j) - <a class="code" href="group__special__functions.html#ga177b023eb871927fcfc90259e9f00382">lgamma</a>(-a * j);</div><div class="line"><a name="l00070"></a><span class="lineno"> 70</span>&#160; ww_max = std::fmax( ww_max, asDouble(ww[k]) );</div><div class="line"><a name="l00071"></a><span class="lineno"> 71</span>&#160; }</div><div class="line"><a name="l00072"></a><span class="lineno"> 72</span>&#160; <span class="keywordflow">for</span> (<span class="keywordtype">int</span> k = 0; k &lt; iterm; k++)</div><div class="line"><a name="l00073"></a><span class="lineno"> 73</span>&#160; sum_ww += exp(ww[k] - ww_max);</div><div class="line"><a name="l00074"></a><span class="lineno"> 74</span>&#160; Float ans = log(sum_ww) + ww_max ;</div><div class="line"><a name="l00075"></a><span class="lineno"> 75</span>&#160;</div><div class="line"><a name="l00076"></a><span class="lineno"> 76</span>&#160; <span class="keywordflow">return</span> ans;</div><div class="line"><a name="l00077"></a><span class="lineno"> 77</span>&#160;}</div><div class="ttc" id="group__special__functions_html_ga177b023eb871927fcfc90259e9f00382"><div class="ttname"><a href="group__special__functions.html#ga177b023eb871927fcfc90259e9f00382">lgamma</a></div><div class="ttdeci">Type lgamma(Type x)</div><div class="ttdoc">Logarithm of gamma function (following R argument convention). </div><div class="ttdef"><b>Definition:</b> <a href="lgamma_8hpp_source.html#l00012">lgamma.hpp:12</a></div></div>
</div><!-- fragment --></div><!-- contents -->
License: <a href="https://gnu.org/licenses/old-licenses/gpl-2.0.txt">GPL v2</a>

0 comments on commit 4ac585c

Please sign in to comment.