Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Unimodality constraint #38

Open
Bio-data-tricks opened this issue Feb 5, 2021 · 5 comments
Open

Unimodality constraint #38

Bio-data-tricks opened this issue Feb 5, 2021 · 5 comments

Comments

@Bio-data-tricks
Copy link

Hello,
I would like to use a Unimodality constraint to follow a kinetic process (time-ordered matrix ). Is it possible to perform this kind of study with your package ?

Thank you

Nicolas

@CCampJr
Copy link
Collaborator

CCampJr commented Feb 5, 2021

Hi Nicolas,

Yes, you can do unimodality with a little bit of coding. In your case, does the mode need to follow an analytical shape or is it just enforcing that the global and local maximum are the same?

For enforcing an analytical modal shape, here's an example: https://pages.nist.gov/pyMCR/usage.html#single-gaussian-distribution-spectral-fitting-with-lmfit

For just enforcing unimodality with no particular shape requirements, I could code up something quickly that you could copy-and-paste until I can push it to the full package. (Unimodality used to be in an older version of the code, but I can't find it -- d'oh!).

Thanks!

@Bio-data-tricks
Copy link
Author

Bio-data-tricks commented Feb 6, 2021 via email

@Bio-data-tricks
Copy link
Author

Bio-data-tricks commented Feb 22, 2021 via email

@CCampJr
Copy link
Collaborator

CCampJr commented Feb 22, 2021

@colapili

Hi Nicolas,

That's a great question and one that I don't think has a singular answer. I have seen it result from preprocessing steps, I've seen it when some of the data are lots of 0's (exactly 0), and it can depend on the type of regressor used.

The first thing I would suspect is that regular/ordinary least-squares (OLS) and non-negative least squares (NNLS) are most likely to cause a noise-like appearance. This is because they will tend to give some weight (regression coefficients) to each and every feature. I would try using Ridge regression (https://pages.nist.gov/pyMCR/usage.html#ridge-regression-from-scikit-learn) with a small weight (alpha) to start -- maybe 1e-5 and continue to lower it by order-of-magnitude until the noise returns. You may even have to raise the alpha to 1 or 10, but I would start small. If Ridge regression is going to converge to OLS, it's usually at 1e-8 to 1e-10 in my experience -- though it doesn't always give the same answer (a case of numerical vs theoretical "should").

Just a side-note: I have found the Ridge SVD solver to be more consistent if using a LOT of data -- maybe your data isn't that big.

Ridge(alpha=1e-5, solver='svd')

Related to that, I can't tell you whether to make the S or the C regressors (or both) Ridge -- just going to need to test it out.

I hope that helps!

Charlie

@CCampJr
Copy link
Collaborator

CCampJr commented Feb 22, 2021

@colapili

Here is an example function that can be used to make a unimodal constraint (see if you can get it going using this example: https://pages.nist.gov/pyMCR/usage.html#single-gaussian-distribution-spectral-fitting-with-lmfit as a template)


import numpy as np

import matplotlib.pyplot as plt
%matplotlib inline

x = np.arange(1001)-500
y = np.exp(-x**2/20**2) + 1.2*np.exp(-(x-100)**2/200**2)
y = np.vstack((y, 2*y))

def unimodal(signal, axis=-1):
    """ Enforces unimodality assuming the peak is upwards"""
    
    assert axis == -1  # Too lazy to implement other currently
    
    if signal.ndim == 1:
        signal = signal[None,:]
        
    loc_max = np.argmax(signal, axis=-1)
    
    # Gradients
    dx_signal_l = np.gradient(signal, axis=axis)  # Left-side of peak
    dx_signal_r = 1*dx_signal_l  # Right-side of peak
    
    for num, lm in enumerate(loc_max):
        dx_signal_l[num,lm:] = 0
        dx_signal_r[num,:lm] = 0
    
        dx_signal_l[num, dx_signal_l[num,:]<0] = 0
        dx_signal_r[num, dx_signal_r[num,:]>0] = 0
    
    return np.cumsum(dx_signal_l + dx_signal_r, axis=-1)
    
    
plt.plot(x,y.T)
plt.plot(x,unimodal(y).T)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants