# Log-scale transform

**URL:** <https://discourse.pymc.io/t/log-scale-transform/12497>\
**Category:** v5\
**Tags:** development, theano, prior, modeling, pytensor\
**Created:** [July 11, 2023, 9:40pm UTC](https://discourse.pymc.io/t/log-scale-transform/12497 "2023-07-11T21:40:17Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![Varun\_Gupta](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/varun_gupta/32/6345_2.png) [@Varun\_Gupta](https://discourse.pymc.io/u/Varun_Gupta)\
**Post date:** [July 11, 2023, 9:40pm UTC](https://discourse.pymc.io/t/log-scale-transform/12497/1 "2023-07-11T21:40:17Z")

</div>

Hi Everyone !

I am implementing a custom likelihood function and doing the computations in the log space.  
There for posterior = prior\*likelihood gets converted to Log posterior= log prior + log likelihood

How do I achieve the log transformation of the priors in pymc? Do I even need to do it explicitly or pymc does it internally?

Thank you !

---

<div class="post-metadata">

**Author:** ![cluhmann](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/cluhmann/32/3083_2.png) [@cluhmann](https://discourse.pymc.io/u/cluhmann)\
**Post date:** [July 12, 2023, 4:39am UTC](https://discourse.pymc.io/t/log-scale-transform/12497/2 "2023-07-12T04:39:18Z")

</div>

Have you checked out the 2 blackbox likelihood notebooks ([the one](https://www.pymc.io/projects/examples/en/latest/case_studies/blackbox_external_likelihood_numpy.html) and [the other](https://www.pymc.io/projects/examples/en/latest/case_studies/blackbox_external_likelihood.html))?

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [July 12, 2023, 7:22am UTC](https://discourse.pymc.io/t/log-scale-transform/12497/3 "2023-07-12T07:22:17Z")

</div>

> [@Varun\_Gupta](#):
>
> How do I achieve the log transformation of the priors in pymc? Do I even need to do it explicitly or pymc does it internally?

You don’t have to worry about the logp of other variables. PyMC will retrieve the logp of each variable (free or observed) and combine it as needed depending on what it is doing.

---

<div class="post-metadata">

**Author:** ![Varun\_Gupta](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/varun_gupta/32/6345_2.png) [@Varun\_Gupta](https://discourse.pymc.io/u/Varun_Gupta)\
**Post date:** [July 12, 2023, 2:36pm UTC](https://discourse.pymc.io/t/log-scale-transform/12497/4 "2023-07-12T14:36:20Z")

</div>

Thanks for referring these, I’ll take a look !

---

<div class="post-metadata">

**Author:** ![Varun\_Gupta](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/varun_gupta/32/6345_2.png) [@Varun\_Gupta](https://discourse.pymc.io/u/Varun_Gupta)\
**Post date:** [July 12, 2023, 2:36pm UTC](https://discourse.pymc.io/t/log-scale-transform/12497/5 "2023-07-12T14:36:33Z")

</div>

Got it ! Thank you !

---

<div class="post-metadata">

**Author:** ![Varun\_Gupta](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/varun_gupta/32/6345_2.png) [@Varun\_Gupta](https://discourse.pymc.io/u/Varun_Gupta)\
**Post date:** [July 13, 2023, 1:28am UTC](https://discourse.pymc.io/t/log-scale-transform/12497/6 "2023-07-13T01:28:02Z")

</div>

Hi @ricardoV94  
Could you also tell me whats the logP function supposed to return?

Log Likelihood computed and summed over all the dataset (this would be one single floating point value OR we could say a scalar) or just likelihood computed for the whole data, without summing (this would be tensor of shape Nx1 where N is the number of records I have in the dataset)

The above calculation is in log space

---

<div class="post-metadata">

**Author:** ![cluhmann](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/cluhmann/32/3083_2.png) [@cluhmann](https://discourse.pymc.io/u/cluhmann)\
**Post date:** [July 13, 2023, 2:47am UTC](https://discourse.pymc.io/t/log-scale-transform/12497/7 "2023-07-13T02:47:28Z")

</div>

@ricardoV94 it seems like it may also depend on whether `pm.Potential` is being use of `pm.CustomDist`? I thought `pm.Potential` required a scalar logp. But can `pm.CustomDist` handle a vector/matrix of logps?

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [July 13, 2023, 3:27am UTC](https://discourse.pymc.io/t/log-scale-transform/12497/8 "2023-07-13T03:27:36Z")

</div>

For potentials it doesn’t ultimately matter because all potentials are summed up in the end, see [here](https://github.com/pymc-devs/pymc/blob/f8581ab84ace71c3a82e3b8eed569c63af84208d/pymc/model.py#L847)

The logp for distributions evaluate at a single point, but can include a batch dimensions. I think it’s most clear when you look at the [logp function for multivariate normal](https://github.com/pymc-devs/pymc/blob/f8581ab84ace71c3a82e3b8eed569c63af84208d/pymc/distributions/multivariate.py#L291). The dimensionality of the distribution, `k`, used in the normalizing constant, is taken to the last dimension of `value`. So you give it an `(n,k)` matrix of values and get back `(n,)` logp numbers. I don’t know of a distribution that doesn’t accept batch dimensions – maybe GPs count?

So does that count as `logp` accepting vector/matrix inputs, or does it just count as broadcasting of a (scalar/vector/matrix, depending on the dimensionality of the support) function over batch dims? I think about it as the latter, but maybe that’s wrong.

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [July 13, 2023, 5:48am UTC](https://discourse.pymc.io/t/log-scale-transform/12497/9 "2023-07-13T05:48:38Z")

</div>

What @jessegrabowski said. Logp should have one entry per independent draws (aka batch dimensions). For univariate rvs, `logp.shape == rv.shape`, for vector rvs `logp.shape==rv.shape[:-1]`, for matrix rvs `logp.shape=rv.shape[:-2]` and so on…

---

<div class="post-metadata">

**Author:** ![Varun\_Gupta](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/varun_gupta/32/6345_2.png) [@Varun\_Gupta](https://discourse.pymc.io/u/Varun_Gupta)\
**Post date:** [July 17, 2023, 10:31pm UTC](https://discourse.pymc.io/t/log-scale-transform/12497/10 "2023-07-17T22:31:14Z")

</div>

if my model output is multivariate, then logp output should be multivariate as well?

One log probability per index of random variable?

Also, one draw one means? One draw while sampling? LogP should spit one value per draw? The shape out Logp O/P depends on random variable under study?

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [July 17, 2023, 10:56pm UTC](https://discourse.pymc.io/t/log-scale-transform/12497/11 "2023-07-17T22:56:21Z")

</div>

Think about the core case. What’s the smallest RV draw you can make? You should have one logp scalar for such a draw.

For a univariate normal you can make a scalar draw with shape `()`, for a multivariate normal a vector draw with shape `(n,)`. Both draws are associated with a scalar logp with shape `()` because you can’t break the draw “apart”. Anything beyond those are “batch dimensions” and you will have corresponding batch dimensions in the logp.

If you a have a univariate normal with batch shapes `(3, 4)`, and thus shape `(3, 4)` as well, then logp will have shape `(3, 4)`.

A multivariate normal with batch shapes `(3, 4)` will have shape `(3, 4, n)` and logp with shape `(3, 4)`. You basically have one entry per batch dimensions.

Batch dimensions are explained in the pymc dimensionality notebook: [Distribution Dimensionality — PyMC 5.6.1 documentation](https://www.pymc.io/projects/docs/en/stable/learn/core_notebooks/dimensionality.html)

---

<div class="post-metadata">

**Author:** ![Varun\_Gupta](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/varun_gupta/32/6345_2.png) [@Varun\_Gupta](https://discourse.pymc.io/u/Varun_Gupta)\
**Post date:** [July 17, 2023, 11:08pm UTC](https://discourse.pymc.io/t/log-scale-transform/12497/12 "2023-07-17T23:08:14Z")

</div>

Thank you @ricardoV94 for being patient.

I have a univariate random variable in my case. According to what you have suggested the draw would have the shape (). Now what I cannot understand is what LogP should return? Lets say I have N data points, should logP return a tensor of shape NX1, where I am returning likelihood calculation for every data point and while sampling PYMC takes care of summing it up OR does LogP return just one single scalar of shape ().

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/5/56b61c26acec21d5876f99db599617af944db36e.png)

Right now my LogP function generates a tensor variable of shape of the observed data (lets say thats N X 1). Currently, logp values are calculated independently for each observation and are not summed while returning from LogP function. Is that the right way to do it? Or do I need to write the summing up part in logP while returning it?

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [July 18, 2023, 6:26am UTC](https://discourse.pymc.io/t/log-scale-transform/12497/13 "2023-07-18T06:26:28Z")

</div>

> [@Varun\_Gupta](#):
>
> Lets say I have N data points, should logP return a tensor of shape NX1,

Almost, should be shape N, not Nx1.

PyMC will sum whenever needed yes.

---

<div class="post-metadata">

**Author:** ![Varun\_Gupta](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/varun_gupta/32/6345_2.png) [@Varun\_Gupta](https://discourse.pymc.io/u/Varun_Gupta)\
**Post date:** [July 18, 2023, 6:33am UTC](https://discourse.pymc.io/t/log-scale-transform/12497/14 "2023-07-18T06:33:40Z")

</div>

Thanks a lot for the clarification!

---

<div class="post-metadata">

**Author:** ![Varun\_Gupta](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/varun_gupta/32/6345_2.png) [@Varun\_Gupta](https://discourse.pymc.io/u/Varun_Gupta)\
**Post date:** [July 18, 2023, 6:25pm UTC](https://discourse.pymc.io/t/log-scale-transform/12497/15 "2023-07-18T18:25:52Z")

</div>

@ricardoV94 I was looking at the example below

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/4/4f2174b227320c6dac88c7ae21d8cacb1d6497f1.png)

Why have they summed while returning from the logp function?

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [July 18, 2023, 6:29pm UTC](https://discourse.pymc.io/t/log-scale-transform/12497/16 "2023-07-18T18:29:22Z")

</div>

That example is unfortunately very outdated.

In general it may also be fine to sum the logp. But some use cases will not be possible if you do this (such as model comparison via elementwise log-likelihood or use in Mixtures, to name a few). So better to do it correctly and not sum unnecessarily.

---

<div class="post-metadata">

**Author:** ![Varun\_Gupta](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/varun_gupta/32/6345_2.png) [@Varun\_Gupta](https://discourse.pymc.io/u/Varun_Gupta)\
**Post date:** [July 18, 2023, 6:39pm UTC](https://discourse.pymc.io/t/log-scale-transform/12497/17 "2023-07-18T18:39:32Z")

</div>

```auto
from pymc.math import where
import pandas as pd
import pymc as pm
from pymc.math import log,exp
from pytensor import tensor as pt
import numpy as np
import nutpie

from numpy import pi
from sklearn.preprocessing import StandardScaler
import pymc.sampling_jax
import numpyro
import jax

numpyro.set_host_device_count(4)

def standardize(x, mu, sigma):
    return (x - mu) / sigma

def standardNormCdf(x):
    return 0.5 + 0.5 * pm.math.erf(x / pm.math.sqrt(2))

def getContributionFromInterval(interval, mu, sigma):
    logT = log(interval)
    a = standardize(logT, mu, sigma)
    return log(1 - standardNormCdf(a))

def computeFailureRate(sigma, mu, t, theta, data):
    logT = log(t)
    a1 = standardize(logT, mu, sigma)
    a2 = 0.5 * log(2 * pi * pow(sigma, 2))
    B = 0.5 * pow(a1, 2)
    C = log(1 - standardNormCdf(a1))
    D = (theta * data).sum(axis=-1).reshape((-1,1))
    failureRate = ((-1) * logT) - a2 - B - C + D
    return failureRate

def getSurvival(time_to_event, mu, sigma, theta, features):
    contributionFromFinalInterval = (-1) * getContributionFromInterval(time_to_event, mu, sigma)
    theta_features_vector = (theta * features).sum(axis=-1)
    theta_features_vector = theta_features_vector.reshape((-1, 1))
    survival = (contributionFromFinalInterval * exp(theta_features_vector))
    # survival = survival.sum(axis=1).reshape((-1, 1))
    return survival

def runModel():
    bcphm_model = pm.Model()
    with bcphm_model:
        data = pd.read_csv('data/precovid_SF.csv')

        data.sort_values(by=['event'], inplace=True)
        data.loc[(data.FirstTimeHomeBuyer.isna()) & data.LoanPurpose.isin(
            ['Refinance', 'CashOutRefi']), 'FirstTimeHomeBuyer'] = 'Not Applicable'
        data['PMI'].fillna(0, inplace=True)
        data = data[data.time_to_event != 0]

        data.dropna(inplace=True)
        data['ClosingDt'] = pd.to_datetime(data['ClosingDt']).dt.year

        data = data.groupby('event', group_keys=False).apply(lambda x: x.sample(frac=0.2))

        events = np.where(data.event.values == "default", 0, np.where(data.event.values == "prepayment", 1, 2)).reshape(
            (-1, 1))

        time_to_event = data.time_to_event.values

        time_to_event_shape = time_to_event.shape
        time_to_event = time_to_event.reshape(time_to_event_shape[0], 1)
        data['isSingleBorrower'] = data['isSingleBorrower'].map({0: 'No', 1: 'Yes'})
        features = data.drop(
            columns=['LoanNumber', 'time_to_default', 'time_to_prepayment', 'State', 'ClosingDt', 'event',
                     'time_to_event'], axis=1)
        categorical = [col for col in features.columns if features[col].dtype == "O"]
        quantitative = set(features.columns) - set(categorical)
        # features = np.repeat(features[:, np.newaxis, :], lifetime.shape[1], axis=1)

        categorical_dummies = pd.get_dummies(features.loc[:, list(categorical)], columns=categorical, drop_first=True)

        sc = StandardScaler()
        standardized_quantitative = pd.DataFrame(sc.fit_transform(features.loc[:, list(quantitative)]),
                                                 columns=list(quantitative))
        model_input = pd.concat(
            [categorical_dummies.reset_index(drop=True), standardized_quantitative.reset_index(drop=True)],
            axis=1)

        model_input.replace({False: 0, True: 1}, inplace=True)

        features_shape = model_input.shape

        features_num = model_input.values
        features_num = pm.MutableData('features_num', features_num)
        events = pm.MutableData('events', events)

        theta_D = pm.Normal('theta_D', mu=0, sigma=100, shape=features_shape[1])
        theta_P = pm.Normal('theta_P', mu=0, sigma=100, shape=features_shape[1])
        mu_D = pm.Normal('mu_D', mu=0, sigma=10)
        mu_P = pm.Normal('mu_P', mu=0, sigma=10)
        sigma_D = pm.Exponential('sigma_D', .01)
        sigma_P = pm.Exponential('sigma_P', .01)

        def logp(time_to_event, mu_P, mu_D, sigma_P,
                 sigma_D, theta_D, theta_P, event,
                 features):
            failureRate = where(
                pt.eq(event, 0),
                computeFailureRate(sigma_D, mu_D, time_to_event, theta_D, features),
                where(pt.eq(event, 1),
                      computeFailureRate(sigma_P, mu_P, time_to_event, theta_P, features),
                      0)
            )
            defaultSurvival = getSurvival(time_to_event, mu_D, sigma_D,
                                          theta_D, features)
            prepaymentSurival = getSurvival(time_to_event, mu_P,
                                            sigma_P, theta_P, features)
            return (failureRate - defaultSurvival - prepaymentSurival).flatten()

       
        likelihood = pm.CustomDist('LL', mu_P, mu_D,
                                   sigma_P, sigma_D, theta_D, theta_P,
                                   events, features_num,
                                   logp=logp,
                                   observed=time_to_event)

    compiled_model = nutpie.compile_pymc_model(bcphm_model)
    trace_pymc = nutpie.sample(compiled_model)
    return trace

```

@ricardoV94 I am having trouble in getting this sampling to converge. Could you please take a look and see any obvious issues? Surprisingly it works on my local machine with vanilla sampler and with smaller fragment of dataset.

Same code with Numpyro on aws sagemaker gives me samples where the variables are stuck and every proposal is rejected.