text:
# Splitters#

Classes that help performing cross-validation.

Our splitters attempt to reduce any potential for data leakage by using grouping by default– and prioritizing grouping over stratficiation or exactly matching the requested train test ratio.

Warning

Due to the grouping operations, the train/test ratios the methods produce will not exactly match the one you requested. For this reason, please get the length of the train/test/valid indices the methods produce.

class DensitySplitter(ds, density_q=None, shuffle=True, random_state=None, sample_frac=1.0, stratification_col=None, center=<function median>, q=(0, 0.25, 0.5, 0.75, 1), sort_by_len=True)[source]#

Splitter that uses the density of the structures to split the data.

For this, we sort structures according to their density and then group the based on the density. You can modify the number of groups using the density_q parameter, those values indicate the quantiles which we use for the grouping.

This ensures that the validation is quite stringent as the different folds will have different densities.

The motivations for doing this are:

• density is often one of the most important descriptors for gas uptake properties.

• there is often is a very large difference in density distribution

between hypothetical and experimental databases.

Initialize the DensitySplitter class.

Parameters:
• ds (AbstractStructureDataset) – A structure dataset. The BaseSplitter only requires the length magic method to be implemented. However, other splitters might require additional methods.

• density_q (Collection[float], optional) – List of quantiles used for quantile binning for the density. Defaults to None. If None, then we use two bins for test/train split, three for validation/train/test split and k for k-fold.

• shuffle (bool) – If True, perform a shuffled split. Defaults to True.

• random_state (Union[int, np.random.RandomState], optional) – Random state for the shuffling. Defaults to None.

• sample_frac (float, optional) – This can be used for downsampling. It will randomly select a subset of indices from all indices before splittings. For instance sample_frac=0.8 will randomly select 80% of the indices before splitting. Defaults to 1.0.

• stratification_col (Union[str, np.typing.ArrayLike], optional) – Data used for stratification. If it is categorical (see mofdscribe.splitters.utils.is_categorical()) then we directly use it for stratification. Otherwise, we use quantile binning. Defaults to None.

• center (callable) – Aggregation function to compute a measure of centrality of all the points in a group such that this can then be used for stratification. This is only used for continuos inputs. For categorical inputs, we always use the mode. Defaults to np.median.

• q (Collection[float], optional) – List of quantiles used for quantile binning. Defaults to (0, 0.25, 0.5, 0.75, 1]. Defaults to [0, 0.25, 0.5, 0.75, 1).

• sort_by_len (bool) – If True, sort the splits by length. (Applies to the train/test/valid and train/test splits). Defaults to True.

class HashSplitter(ds, hash_type='undecorated_scaffold_hash', shuffle=True, random_state=None, sample_frac=1.0, stratification_col=None, center=<function median>, q=(0, 0.25, 0.5, 0.75, 1), sort_by_len=True)[source]#

Splitter that uses Weisfeiller-Lehman graph hashes [WL] to split the data in more stringent ways.

Note that the hashes we use do not allow for a meaningful measure of similarity. That is, there is no way to measure the distance between two strings. The only meaningful measure is if they are identical or not.

Note

Weisfeiller-Lehman graph hashes do not give a guarantee for graph-isomorphism. That is, there might be identical hashes that do not correspond to isomorphic graphs.

Note

There are certain graphs that a Weisfeiller-Lehman test cannot distinguish [Bouritsas].

Note

We speak about Weisfeiller-Lehman hashes as they are the defaults for the mofdscribe datasets. However, you can also overwrite this method with a custom hashing function.

Initialize a HashSplitter.

Parameters:
• ds (AbstractStructureDataset) – A structure dataset. The BaseSplitter only requires the length magic method to be implemented. However, other splitters might require additional methods.

• hash_type (str) – Hash type to use. Must be one of the following: * undecorated_scaffold_hash * decorated_graph_hash * decorated_scaffold_hash * undecorated_graph_hash Defaults to “undecorated_scaffold_hash”.

• shuffle (bool) – If True, perform a shuffled split. Defaults to True.

• random_state (Union[int, np.random.RandomState], optional) – Random state for the shuffling. Defaults to None.

• sample_frac (float, optional) – This can be used for downsampling. It will randomly select a subset of indices from all indices before splittings. For instance sample_frac=0.8 will randomly select 80% of the indices before splitting. Defaults to 1.0.

• stratification_col (Union[str, np.typing.ArrayLike], optional) – Data used for stratification. If it is categorical (see mofdscribe.splitters.utils.is_categorical()) then we directly use it for stratification. Otherwise, we use quantile binning. Defaults to None.

• center (callable, optional) – Aggregation function to compute a measure of centrality of all the points in a group such that this can then be used for stratification. This is only used for continuos inputs. For categorical inputs, we always use the mode. Defaults to np.median.

• q (Collection[float], optional) – List of quantiles used for quantile binning. Defaults to (0, 0.25, 0.5, 0.75, 1).

• sort_by_len (bool) – If True, sort the splits by length. (Applies to the train/test/valid and train/test splits). Defaults to True.

class TimeSplitter(ds, year_q=None, shuffle=True, random_state=None, sample_frac=1.0, stratification_col=None, center=<function median>, q=(0, 0.25, 0.5, 0.75, 1), sort_by_len=True)[source]#

This splitter sorts structures according to their publication date.

That is, the training set will contain structures that are “older” (have been discovered earlier) than the ones in the test set. This can mimick real-life model development conditions [MoleculeNet].

It has for instance also be used with ICSD data in [Palizhati] and been the focus of [Sheridan].

Initialize the TimeSplitter class.

Parameters:
• ds (AbstractStructureDataset) – A structure dataset. The BaseSplitter only requires the length magic method to be implemented. However, other splitters might require additional methods.

• year_q (Collection[float]) – List of quantiles used for quantile binning on the years. Defaults to None. If None, then we use two bins for test/train split, three for validation/train/test split and k for k-fold.

• shuffle (bool) – If True, perform a shuffled split. Defaults to True.

• random_state (Union[int, np.random.RandomState], optional) – Random state for the shuffling. Defaults to None.

• sample_frac (float, optional) – This can be used for downsampling. It will randomly select a subset of indices from all indices before splittings. For instance sample_frac=0.8 will randomly select 80% of the indices before splitting. Defaults to 1.0.

• stratification_col (Union[str, np.typing.ArrayLike], optional) – Data used for stratification. If it is categorical (see mofdscribe.splitters.utils.is_categorical()) then we directly use it for stratification. Otherwise, we use quantile binning. Defaults to None.

• center (callable) – Aggregation function to compute a measure of centrality of all the points in a group such that this can then be used for stratification. This is only used for continuos inputs. For categorical inputs, we always use the mode. Defaults to np.median.

• q (Collection[float], optional) – List of quantiles used for quantile binning. Defaults to (0, 0.25, 0.5, 0.75, 1).

• sort_by_len (bool) – If True, sort the splits by length. (Applies to the train/test/valid and train/test splits). Defaults to True.

class BaseSplitter(ds, shuffle=True, random_state=None, sample_frac=1.0, stratification_col=None, center=<function median>, q=(0, 0.25, 0.5, 0.75, 1), sort_by_len=True)[source]#

A BaseSplitter implements the basic logic for dataset partition as well as k-fold cross-validation.

Methods that inherit from this class typically implement the

• code:

_get_stratification_col: Should return an ArrayLike object of floats, categories, or ints. If it is categorical data, the BaseSplitter will handle the discretization.

• code:

_get_groups: Should return an ArrayLike object of categories (integers or strings)

methods. Internally, the BaseSplitter uses those to group and/or stratify the splits.

Initialize a BaseSplitter.

Parameters:
• ds (AbstractStructureDataset) – A structure dataset. The BaseSplitter only requires the length magic method to be implemented. However, other splitters might require additional methods.

• shuffle (bool) – If True, perform a shuffled split. Defaults to True.

• random_state (Optional[Union[int, np.random.RandomState]], optional) – Random state for the shuffling. Defaults to None.

• sample_frac (Optional[float], optional) – This can be used for downsampling. It will randomly select a subset of indices from all indices before splittings. For instance sample_frac=0.8 will randomly select 80% of the indices before splitting. Defaults to 1.0.

• stratification_col (Union[str, np.typing.ArrayLike], optional) – Data used for stratification. If it is categorical (see mofdscribe.splitters.utils.is_categorical()) then we directly use it for stratification. Otherwise, we use quantile binning. Defaults to None.

• center (callable) – Aggregation function to compute a measure of centrality of all the points in a group such that this can then be used for stratification. This is only used for continuos inputs. For categorical inputs, we always use the mode. Defaults to np.median.

• q (Collection[float], optional) – List of quantiles used for quantile binning. Defaults to (0, 0.25, 0.5, 0.75, 1).

• sort_by_len (bool) – If True, sort the splits by length. (Applies to the train/test/valid and train/test splits). Defaults to True.

train_test_split(frac_train=0.7)[source]#

Perform a train/test partition.

Parameters:

frac_train (float) – Fraction of the data to use for the training set. Defaults to 0.7.

Returns:

Train indices, test indices

Return type:

Tuple[Collection[int], Collection[int]]

train_valid_test_split(frac_train=0.7, frac_valid=0.1)[source]#

Perform a train/valid/test partition.

Parameters:
• frac_train (float) – Fraction of data to use for the training set. Defaults to 0.7.

• frac_valid (float) – Fraction of data to use for the validation set. Defaults to 0.1.

Returns:

Training, validation, test set.

Return type:

Tuple[Collection[int], Collection[int], Collection[int]]

k_fold(k=5)[source]#

Peform k-fold crossvalidation.

Parameters:

k (int) – Number of folds. Defaults to 5.

Yields:

Tuple[Collection[int], Collection[int]] – Train indices, test indices.

Return type:

Tuple[Collection[int], Collection[int]]

class KennardStoneSplitter(ds, feature_names, shuffle=True, random_state=None, sample_frac=1.0, scale=True, centrality_measure='mean', metric='euclidean', ascending=False)[source]#

Run the Kennard-Stone sampling algorithm [KennardStone].

The algorithm selects samples with uniform converage. The initial samples are biased towards the boundaries of the dataset. Hence, it might be biased by outliers.

This algorithm ensures a flat coverage of the dataset. It is also known as CADEX algorithm and has been later refined in the DUPLEX algorithm [Snee].

Warning

This splitter can be slow for large datasets as it requires us to perform distance matrices N times for a dataset with N structures.

Warning

Stratification is not supported for this splitter.

Warning

I couldn’t find a good reference for the k-fold version of this algorihm.

Construct a KennardStoneSplitter.

Parameters:
• ds (AbstractStructureDataset) – A structure dataset. The BaseSplitter only requires the length magic method to be implemented. However, other splitters might require additional methods.

• feature_names (List[str]) – Names of features to consider.

• shuffle (bool) – If True, perform a shuffled split. Defaults to True.

• random_state (Union[int, np.random.RandomState], optional) – Random state for the shuffling. Defaults to None.

• sample_frac (float, optional) – This can be used for downsampling. It will randomly select a subset of indices from all indices before splittings. For instance sample_frac=0.8 will randomly select 80% of the indices before splitting. Defaults to 1.0.

• scale (bool) – If True, apply z-score normalization prior to running the sampling. Defaults to True.

• centrality_measure (str) – The first sample is selected to be maximally distanct from this value. It can be one of “mean”, “median”, “random”. In case of “random” we simply select a random point. In the case of “mean” and “median” the initial point is maximally distanct from the mean and median of the feature matrix, respectively. Defaults to “mean”.

• metric (Union[Callable, str]) – The distance metric to use. If a string, the distance function can be ‘braycurtis’, ‘canberra’, ‘chebyshev’, ‘cityblock’, ‘correlation’, ‘cosine’, ‘dice’, ‘euclidean’, ‘hamming’, ‘jaccard’, ‘jensenshannon’, ‘kulsinski’, ‘kulczynski1’, ‘mahalanobis’, ‘matching’, ‘minkowski’, ‘rogerstanimoto’, ‘russellrao’, ‘seuclidean’, ‘sokalmichener’, ‘sokalsneath’, ‘sqeuclidean’, ‘yule’. Defaults to “euclidean”.

• ascending (bool) – If True, sort samples in asceding distance to the center. That is, the first samples (maximally distant to center) would be sampled last. Defaults to False.

get_sorted_indices(ds)[source]#

Return a list of indices, sorted by similarity using the Kennard-Stone algorithm.

The first sample will be maximally distant from the center.

Parameters:

ds (AbstractStructureDataset) – A mofdscribe AbstractStructureDataset

Returns:

Sorted indices.

Return type:

Collection[int]

train_test_split(frac_train=0.7)[source]#

Perform a train/test partition.

Parameters:

frac_train (float) – Fraction of the data to use for the training set. Defaults to 0.7.

Returns:

Train indices, test indices

Return type:

Tuple[Collection[int], Collection[int]]

train_valid_test_split(frac_train=0.7, frac_valid=0.1)[source]#

Perform a train/valid/test partition.

Parameters:
• frac_train (float) – Fraction of data to use for the training set. Defaults to 0.7.

• frac_valid (float) – Fraction of data to use for the validation set. Defaults to 0.1.

Returns:

Training, validation, test set.

Return type:

Tuple[Collection[int], Collection[int], Collection[int]]

k_fold(k=5)[source]#

Peform k-fold crossvalidation.

Parameters:

k (int) – Number of folds. Defaults to 5.

Yields:

Tuple[Collection[int], Collection[int]] – Train indices, test indices.

Return type:

Tuple[Collection[int], Collection[int]]

class ClusterSplitter(ds, feature_names, shuffle=True, random_state=None, sample_frac=1.0, stratification_col=None, center=<function median>, q=(0, 0.25, 0.5, 0.75, 1), sort_by_len=False, scaled=True, n_pca_components='mle', n_clusters=4, pca_kwargs=None, kmeans_kwargs=None)[source]#

Split the data into clusters and use the clusters as groups.

The approach has been proposed on Kaggle. In principle, we perform the following steps:

1. Scale the data (optional).

2. Perform PCA for de-correlation.

3. Perform k-means clustering.

Construct a ClusterSplitter.

Parameters:
• ds (AbstractStructureDataset) – A structure dataset. The BaseSplitter only requires the length magic method to be implemented. However, other splitters might require additional methods.

• feature_names (List[str]) –

Names of features to consider.

shuffle (bool): If True, perform a shuffled split.

Defaults to True.

• shuffle (bool) – If True, perform a shuffled split. Defaults to True.

• random_state (int) – Random state for the shuffling. Defaults to None.

• sample_frac (float, optional) – This can be used for downsampling. It will randomly select a subset of indices from all indices before splittings. For instance sample_frac=0.8 will randomly select 80% of the indices before splitting. Defaults to 1.0.

• stratification_col (Union[str, np.typing.ArrayLike], optional) – Data used for stratification. If it is categorical (see mofdscribe.splitters.utils.is_categorical()) then we directly use it for stratification. Otherwise, we use quantile binning. Defaults to None.

• center (callable) – Aggregation function to compute a measure of centrality of all the points in a group such that this can then be used for stratification. This is only used for continuos inputs. For categorical inputs, we always use the mode. Defaults to np.median.

• q (Collection[float], optional) – List of quantiles used for quantile binning. Defaults to (0, 0.25, 0.5, 0.75, 1]. Defaults to [0, 0.25, 0.5, 0.75, 1).

• sort_by_len (bool) – If True, sort the splits by length. (Applies to the train/test/valid and train/test splits). Defaults to True.

• scaled (bool) – If True, scale the data before clustering. Defaults to True.

• n_pca_components (Union[int, str]) – Number of components to use for PCA. If “mle”, use the number of components that maximizes the variance. Defaults to “mle”.

• n_clusters (int) – Number of clusters to use. Defaults to 4.

• random_state – Random seed. Defaults to 42.

• pca_kwargs (Dict[str, Any]) – Keyword arguments to pass to PCA. Defaults to None.

• kmeans_kwargs (Dict[str, Any]) – Keyword arguments to pass to k-means. Defaults to None.

class LOCOCV(ds, feature_names, shuffle=True, random_state=None, sample_frac=1.0, scaled=True, n_pca_components='mle', pca_kwargs=None, kmeans_kwargs=None)[source]#

Leave-one-cluster-out cross-validation.

The general idea has been discussed before, e.g. in [Kramer]. Perhaps more widely used in the materials community is [Meredig]. Here, we perform PCA, followed by k-means clustering.

• Where k = 2 for a train/test split

• Where k = 3 for a train/valid/test split

• Where k = k for k-fold crossvalidation

By default, we will sort outputs such that the cluster sizes are train >= test >= valid.

Construct a LOCOCV.

Parameters:
• ds (AbstractStructureDataset) – A structure dataset. The BaseSplitter only requires the length magic method to be implemented. However, other splitters might require additional methods.

• feature_names (List[str]) – Names of features to consider.

• shuffle (bool) – If True, perform a shuffled split. Defaults to True.

• random_state (int) – Random state for the shuffling. Defaults to None.

• sample_frac (float, optional) – This can be used for downsampling. It will randomly select a subset of indices from all indices before splittings. For instance sample_frac=0.8 will randomly select 80% of the indices before splitting. Defaults to 1.0.

• scaled (bool) – If True, scale the data before clustering. Defaults to True.

• n_pca_components (int) – Number of components to use for PCA. If “mle”, use the number of components that maximizes the variance. Defaults to “mle”.

• random_state – Random seed. Defaults to 42.

• pca_kwargs (Dict[str, Any], optional) – Additional keyword arguments for sklearn’s sklearn.decomposition.PCA. Defaults to None.

• kmeans_kwargs (Dict[str, Any], optional) – Additional keyword arguments for sklearn’s sklearn.clustering.KMeans. Defaults to None.

train_test_split()[source]#

Perform a train/test partition.

Returns:

Train indices, test indices

Return type:

Tuple[Collection[int], Collection[int]]

train_valid_test_split()[source]#

Perform a train/valid/test partition.

Returns:

Training, validation, test set.

Return type:

Tuple[Collection[int], Collection[int], Collection[int]]

k_fold(k)[source]#

Peform k-fold crossvalidation.

Parameters:

k (int) – Number of folds. Defaults to 5.

Yields:

Iterator[Tuple[Collection[int], Collection[int]]] – Train indices, test indices.

Return type:

Tuple[Collection[int], Collection[int]]

class ClusterStratifiedSplitter(ds, feature_names, shuffle=True, random_state=None, sample_frac=1.0, scaled=True, n_pca_components='mle', n_clusters=4, pca_kwargs=None, kmeans_kwargs=None)[source]#

Split the data into clusters and stratify on those clusters

The approach has been proposed on Kaggle. In principle, we perform the following steps:

1. Scale the data (optional).

2. Perform PCA for de-correlation.

3. Perform k-means clustering.

Construct a ClusterStratifiedSplitter.

Parameters:
• ds (AbstractStructureDataset) – A structure dataset. The BaseSplitter only requires the length magic method to be implemented. However, other splitters might require additional methods.

• feature_names (List[str]) – Names of features to consider.

• shuffle (bool) – If True, perform a shuffled split. Defaults to True.

• random_state (int) – Random state for the shuffling. Defaults to None.

• sample_frac (float, optional) – This can be used for downsampling. It will randomly select a subset of indices from all indices before splittings. For instance sample_frac=0.8 will randomly select 80% of the indices before splitting. Defaults to 1.0.

• scaled (bool) – If True, scale the data before clustering. Defaults to True.

• n_pca_components (int) – Number of components to use for PCA. If “mle”, use the number of components that maximizes the variance. Defaults to “mle”.

• n_clusters (int) – Number of clusters to use. Defaults to 4.

• random_state – Random seed. Defaults to 42.

• pca_kwargs (Dict[str, Any]) – Keyword arguments to pass to PCA. Defaults to None.

• kmeans_kwargs (Dict[str, Any]) – Keyword arguments to pass to k-means. Defaults to None.

Helper functions for the splitters.

Some of these methods might also be useful for constructing nested cross-validation loops.

kennard_stone_sampling(X, scale=True, centrality_measure='mean', metric='euclidean')[source]#

Run the Kennard-Stone sampling algorithm [KennardStone].

The algorithm selects samples with uniform converage. The initial samples are biased towards the boundaries of the dataset.

Note

You also might this algorithm useful for creating a “diverse” sample of points, e.g., to initalize an active learnign loop.

Warning

This algorithm has a high computational complexity. It is not recommended for large datasets.

Parameters:
• X (ArrayLike) – Input feature matrix.

• scale (bool) – If True, apply z-score normalization prior to running the sampling. Defaults to True.

• centrality_measure (str) – The first sample is selected to be maximally distanct from this value. It can be one of “mean”, “median”, “random”. In case of “random” we simply select a random point. In the case of “mean” and “median” the initial point is maximally distanct from the mean and median of the feature matrix, respectively. Defaults to “mean”.

• metric (Union[Callable, str]) –

The distance metric to use. If a string, the distance function can be ‘braycurtis’, ‘canberra’, ‘chebyshev’, ‘cityblock’, ‘correlation’, ‘cosine’, ‘dice’, ‘euclidean’, ‘hamming’, ‘jaccard’, ‘jensenshannon’, ‘kulsinski’, ‘kulczynski1’, ‘mahalanobis’, ‘matching’, ‘minkowski’, ‘rogerstanimoto’, ‘russellrao’, ‘seuclidean’, ‘sokalmichener’, ‘sokalsneath’, ‘sqeuclidean’, ‘yule’. See scipy.spatial.distance.cdist().

Defaults to “euclidean”.

Raises:

ValueError – If non-implemented centrality measure is used.

Returns:

indices sorted by their Max-Min distance.

Return type:

List[int]

pca_kmeans(X, scaled, n_pca_components, n_clusters, random_state=None, pca_kwargs=None, kmeans_kwargs=None)[source]#

Run principal component analysis (PCA) followed by K-means clustering on the data.

Uses sklearn’s implementation of PCA, and k-means.

Parameters:
• X (np.ndarray) – Input data

• scaled (bool) – If True, use standard scaling for clustering

• n_pca_components (Union[int, str]) – number of principal components to keep

• n_clusters (int) – number of clusters

• random_state (Optional[Union[int, np.random.RandomState]], optional) – Random state for sklearn. Defaults to None.

• pca_kwargs (Dict[str, Any], optional) – Additional keyword arguments for sklearn’s sklearn.decomposition.PCA. Defaults to None.

• kmeans_kwargs (Dict[str, Any], optional) – Additional keyword arguments for sklearn’s sklearn.clustering.KMeans. Defaults to None.

Returns:

Cluster indices.

Return type:

np.ndarray

is_categorical(x)[source]#

Return true if x is categorial or composed of integers.

Return type:

bool

stratified_train_test_partition(idxs, stratification_col, train_size, valid_size, test_size, shuffle=True, random_state=None, q=(0, 0.25, 0.5, 0.75, 1))[source]#

Perform a stratified train/test split.

Parameters:
• idxs (Sequence[int]) – Indices of points to split

• stratification_col (np.typing.ArrayLike) – Data used for stratification. If it is categorical (see mofdscribe.splitters.utils.is_categorical()) then we directly use it for stratification. Otherwise, we use quantile binning.

• train_size (float) – Size of the training set as fraction.

• valid_size (float) – Size of the validation set as fraction.

• test_size (float) – Size of the test set as fraction.

• shuffle (bool) – If True, perform a shuffled split. Defaults to True.

• random_state (Union[int, np.random.RandomState], optional) – Random state for the suffler. Defaults to None.

• q (Sequence[float], optional) – List of quantiles used for quantile binning. Defaults to (0, 0.25, 0.5, 0.75, 1).

Returns:

Train, validation, test indices.

Return type:

Tuple[np.array, np.array, np.array]

grouped_stratified_train_test_partition(stratification_col, group_col, train_size, valid_size, test_size, shuffle=True, random_state=None, q=(0, 0.25, 0.5, 0.75, 1), center=<function median>)[source]#

Return grouped stratified train-test partition.

First, we compute the most common stratification category / centrality measure of the stratification column for every group. Then, we perform a stratified train/test partition on the groups. We then “expand” by concatenating the indices belonging to each group.

Warning

Note that this won’t work well if the number of groups and datapoints is small. It will also cause issues if the number of datapoints in the groups is very imbalanced.

Parameters:
• stratification_col (np.typing.ArrayLike) – Data used for stratification. If it is categorical (see mofdscribe.splitters.utils.is_categorical()) then we directly use it for stratification. Otherwise, we use quantile binning.

• group_col (np.typing.ArrayLike) – Data used for grouping.

• train_size (float) – Size of the training set as fraction.

• valid_size (float) – Size of the validation set as fraction.

• test_size (float) – Size of the test set as fraction.

• shuffle (bool) – If True, perform a shuffled split. Defaults to True.

• random_state (Union[int, np.random.RandomState], optional) – Random state for the suffler. Defaults to None.

• q (Sequence[float], optional) – List of quantiles used for quantile binning. Defaults to [0, 0.25, 0.5, 0.75, 1].

• center (Callable) – Aggregation function to compute a measure of centrality of all the points in a group such that this can then be used for stratification. This is only used for continuos inputs. For categorical inputs, we always use the mode.

Returns:

Train, validation, test indices.

Return type:

Tuple[np.array, np.array, np.array]

get_train_valid_test_sizes(size, train_size, valid_size, test_size)[source]#

Compute the number of points in every split.

Return type:

Tuple[int, int, int]

grouped_train_valid_test_partition(groups, train_size, valid_size, test_size, shuffle=True, random_state=None)[source]#

Perform a grouped train/test split without stratification.

Parameters:
• groups (np.typing.ArrayLike) – Data used for grouping.

• train_size (float) – Size of the training set as fraction.

• valid_size (float) – Size of the validation set as fraction.

• test_size (float) – Size of the test set as fraction.

• shuffle (bool) – If True, perform a shuffled split. Defaults to True.

• random_state (Union[int, np.random.RandomState], optional) – Random state for the suffler. Defaults to None.

Returns:

Train, validation, test indices.

Return type:

Tuple[np.array, Optional[np.array], np.array]

quantile_binning(values, q)[source]#

Use pandas.qcut() to bin the values based on quantiles.

Return type:

array

check_fraction(train_fraction, valid_fraction, test_fraction)[source]#

Check that the fractions are all between 0 and 1 and that they sum up to 1.

Return type:

None

no_group_warn(groups)[source]#

Raise warning if groups is None.

Return type:

None


text:
# Theta function identity

Let us consider the Theta function $$\theta(\tau) = \sum_{n \in \mathbb{Z}} e^{\pi i n^2 \tau} \text{ for } \mathrm{Im}(\tau)>0.$$ Then it is rather easy to see that $$\theta(\tau+2)=\theta(\tau)$$ and $$\theta(\tau)+\theta(\tau+1)=2 \cdot \theta(4 \tau).$$ Further, one can prove that $$\theta(-1/ \tau) = \sqrt{\frac{\tau}{i}} \cdot \theta(\tau).$$ All together, this yields $$\theta(1- \frac{1}{\tau}) = \theta(- \frac{1}{\tau}) + \theta(1- \frac{1}{\tau}) - \theta(- \frac{1}{\tau}) = 2 \cdot \theta(- \frac{4}{\tau}) - \theta(- \frac{1}{\tau}) = \sqrt{\frac{\tau}{i}} \cdot (\theta(\frac{\tau}{4}) - \theta(\tau)).$$ Here $$\sqrt{}$$ always denotes the principal branch of the square root.

Next, my textbook claims that $$f(\tau) := \theta^4(\tau)-\theta^4(\tau+1)+\tau^{-2} \cdot \theta^4(1- \frac{1}{\tau})$$ fulfils both $$f(\tau+1) = - f(\tau)$$ and $$f(- \frac{1}{\tau}) = - \tau^{2} \cdot f(\tau)$$ (and uses this to deduce that actually $$f = 0$$). It should be possible to derive these using the above $$\theta$$-identities, but I am completely stuck here. Any help is appreciated!

Let us compute explicitly: \begin{aligned} f(\tau)+f(\tau+1) &= \theta^4(\tau)-\theta^4(\tau+1) +\tau^{-2} \cdot \theta^4\left(1- \frac{1}{\tau}\right) \\ &\qquad\qquad + \theta^4(\tau+1)-\theta^4(\tau+2) +(\tau+1)^{-2} \cdot \theta^4\left(1- \frac{1}{\tau+1}\right) \\ &= \tau^{-2} \cdot \theta^4\left(1- \frac{1}{\tau}\right) +(\tau+1)^{-2} \cdot \theta^4\left(1- \frac{1}{\tau+1}\right) \\ &= \tau^{-2} \cdot\left[\ \theta^4\left(\frac{\tau-1}{\tau}\right) +\left(\frac\tau{\tau+1}\right)^2 \cdot \theta^4\left(\frac{\tau}{\tau+1}\right) \ \right] \\ &= \tau^{-2} \cdot\left[\ \theta^4\left(\frac{\tau-1}{\tau}-2\right) -\theta^4\left(-\frac 1{\qquad\frac{\tau}{\tau+1}\qquad}\right) \ \right] \\ &=0\ , \\[3mm] f\left(-\frac1\tau\right) &= \theta^4\left(-\frac1\tau\right) - \theta^4\left(-\frac1\tau+1\right) + \left(-\frac1\tau\right)^{-2}\theta^4\left(1-\frac1{-1/\tau}\right) \\ &= \theta^4\left(-\frac1\tau\right) - \theta^4\left(1-\frac1\tau\right) + \tau^2\theta^4(\tau+1) \\ &= -\tau^2\theta^4(\tau) +\tau^2\theta^4(\tau+1) - \tau^2\cdot\tau^{-2} \theta^4\left(1-\frac1\tau\right) \\ &= -\tau^2\;f(\tau)\ . \end{aligned} (The power of $$\tau$$ differs in the last relation.)
• Yes, the power of $\tau$ in the last relation should be $2$. I edited it. Thank you! Sep 3, 2019 at 12:10

format: plain


text:
# Properties

 Label 117600dd Number of curves $4$ Conductor $117600$ CM no Rank $1$ Graph

# Related objects

Show commands for: SageMath
sage: E = EllipticCurve("dd1")

sage: E.isogeny_class()

## Elliptic curves in class 117600dd

sage: E.isogeny_class().curves

LMFDB label Cremona label Weierstrass coefficients j-invariant Discriminant Torsion structure Modular degree Faltings height Optimality
117600.hz3 117600dd1 $$[0, 1, 0, -86158, -9607312]$$ $$601211584/11025$$ $$1297080225000000$$ $$[2, 2]$$ $$884736$$ $$1.6947$$ $$\Gamma_0(N)$$-optimal
117600.hz4 117600dd2 $$[0, 1, 0, -408, -27786312]$$ $$-8/354375$$ $$-333534915000000000$$ $$[2]$$ $$1769472$$ $$2.0413$$
117600.hz2 117600dd3 $$[0, 1, 0, -178033, 14372063]$$ $$82881856/36015$$ $$271176239040000000$$ $$[2]$$ $$1769472$$ $$2.0413$$
117600.hz1 117600dd4 $$[0, 1, 0, -1372408, -619289812]$$ $$303735479048/105$$ $$98825160000000$$ $$[2]$$ $$1769472$$ $$2.0413$$

## Rank

sage: E.rank()

The elliptic curves in class 117600dd have rank $$1$$.

## Complex multiplication

The elliptic curves in class 117600dd do not have complex multiplication.

## Modular form 117600.2.a.dd

sage: E.q_eigenform(10)

$$q + q^{3} + q^{9} + 4q^{11} + 6q^{13} - 6q^{17} - 4q^{19} + O(q^{20})$$

## Isogeny matrix

sage: E.isogeny_class().matrix()

The $$i,j$$ entry is the smallest degree of a cyclic isogeny between the $$i$$-th and $$j$$-th curve in the isogeny class, in the Cremona numbering.

$$\left(\begin{array}{rrrr} 1 & 2 & 2 & 2 \\ 2 & 1 & 4 & 4 \\ 2 & 4 & 1 & 4 \\ 2 & 4 & 4 & 1 \end{array}\right)$$

## Isogeny graph

sage: E.isogeny_graph().plot(edge_labels=True)

The vertices are labelled with Cremona labels.


text:
# Find the general solution for the differential equation

So in my previous math class I spotted on my book an exercise that I couldnt solve. We had to find the general solution for the differential equation. This was the exercise: 4y'' - 4y' + y = ex/2√(1-x2)
Can anyone tell me how to solve this step by step?

fresh_42
Mentor
You first solve the homogeneous version.

HallsofIvy
Homework Helper
The "associated homogeneous equation" is -4y'+ y= 0 so y'= -y/4. We can write that as dy/y= -dx/4. Integrating, ln|y|= -x/4+ C. Taking the exponential of both sides $$y= C' e^{-x/4}$$ where $$|C'|= e^C$$ and we can remove the absolute value by allowing C' to be either positive or negative.
To solve the entire equation use "variation of parameters": Look for a solution of the form $$y(x)= u(x)e^{-x/4}$$. Then $$y'= u'e^{-x/4}- (u/4)e^{-x/4}= (u'- u/4)e^{-x/4}$$ so the equation becomes $$(4u'- u)e^{-x/4}+ ue^{-x/4}= 4u'e^{-x/4}= e^{x/2}\sqrt{1- x^2}$$ so that $$u'= \frac{1}{4}e^{3x/4}(1- x^2)^{1/2}$$. Can you integrate that?

Delta2
Mark44
Mentor
We had to find the general solution for the differential equation. This was the exercise: 4y'' - 4y' + y = ex/2√(1-x2)

The "associated homogeneous equation" is -4y'+ y= 0
No. You apparently missed the second-order term 4y''.

Can anyone tell me how to solve this step by step?
That's not what we do here. Per our rules (https://www.physicsforums.com/threads/physics-forums-global-guidelines.414380/, under Homework Guidelines), you must first make an attempt. Also, homework-type questions should be posted in one of the sections under Homework & Coursework.

HallsofIvy


text:
# Chapter 11 - Counting Methods and Probability Theory - 11.2 Permutations - Exercise Set 11.2 - Page 701: 18

980,179,200

#### Work Step by Step

$\frac{17!}{9!}$ = $\frac{17*16*15*14*13*12*11*10*9!}{9!}$ 9! is common in the numerator and the denominator of the given expression, so cancel 9!. Therefore, evaluate $17*16*15*14*13*12*11*10 = 980,179,200$

After you claim an answer you’ll have 24 hours to send in a draft. An editor will review the submission and either publish your submission or provide feedback.


text:
## Trivial, but useful: sequences with defined mean/s.d.

July 31, 2013

O.k., the following post may be (mathematically) trivial, but could be somewhat useful for people that do simulations/testing of statistical methods.
Let’s say we want to test the dependence of p-values derived from a t-test to a) the ratio of means between two groups, b) the standard deviation or c) the sample size(s) of the two groups. For this setup we would need to i.e. generate two groups with defined $\mu, \sigma$ and $n$.
Often encountered in simulations is that groups are generated with rnorm and then plugged into the simulation. However (and evidently), it is clear that sampling from a normal distribution does not deliver a vector with exactly defined statistical properties (although the “law of large numbers” states that with enough large sample size it converges to that…).
For example,

> x <- rnorm(1000, 5, 2)
> mean(x)
[1] 4.998388
> sd(x)
[1] 2.032262

shows what I meant above ($\mu_x \neq 5, \sigma_x \neq 2$).

Luckily, we can create vectors with exact mean and s.d. by a “scaled-and-shifted z-transformation” of an input vector $X$:

$Z = \frac{X - \mu_X}{\sigma_X} \cdot \mathbf{sd} + \mathbf{mean}$

where sd is the desired standard deviation and mean the desired mean of the output vector Z.

The code is simple enough:

statVec <- function(x, mean, sd)
{
X <- x
MEAN <- mean
SD <- sd
Z <- (((X - mean(X, na.rm = TRUE))/sd(X, na.rm = TRUE))) * SD
MEAN + Z
}

So, using this on the rnorm-generated vector x from above:

> z <- statVec(x, 5, 2)
> mean(z)
[1] 5
> sd(z)
[1] 2

we have created a vector with exact statistical properties, which is also normally distributed since multiplication and addition of a normal distribution preserves normality.

Cheers, Andrej

## wapply: A faster (but less functional) ‘rollapply’ for vector setups

April 23, 2013

For some cryptic reason I needed a function that calculates function values on sliding windows of a vector. Googling around soon brought me to ‘rollapply’, which when I tested it seems to be a very versatile function. However, I wanted to code my own version just for vector purposes in the hope that it may be somewhat faster.
This is what turned out (wapply for “window apply”):

wapply <- function(x, width, by = NULL, FUN = NULL, ...)
{
FUN <- match.fun(FUN)
if (is.null(by)) by <- width

lenX <- length(x)
SEQ1 <- seq(1, lenX - width + 1, by = by)
SEQ2 <- lapply(SEQ1, function(x) x:(x + width - 1))

OUT <- lapply(SEQ2, function(a) FUN(x[a], ...))
OUT <- base:::simplify2array(OUT, higher = TRUE)
return(OUT)
}

It is much more restricted than ‘rollapply’ (no padding, left/center/right adjustment etc).
But interestingly, for some setups it is very much faster:

library(zoo)
x <- 1:200000

large window, small slides:

> system.time(RES1 <- rollapply(x, width = 1000, by = 50, FUN = fun))
User      System verstrichen
3.71        0.00        3.84
> system.time(RES2 <- wapply(x, width = 1000, by = 50, FUN = fun))
User      System verstrichen
1.89        0.00        1.92
> all.equal(RES1, RES2)
[1] TRUE

small window, small slides:

> system.time(RES1 <- rollapply(x, width = 50, by = 50, FUN = fun))
User      System verstrichen
2.59        0.00        2.67
> system.time(RES2 <- wapply(x, width = 50, by = 50, FUN = fun))
User      System verstrichen
0.86        0.00        0.89
> all.equal(RES1, RES2)
[1] TRUE

small window, large slides:

> system.time(RES1 <- rollapply(x, width = 50, by = 1000, FUN = fun))
User      System verstrichen
1.68        0.00        1.77
> system.time(RES2 <- wapply(x, width = 50, by = 1000, FUN = fun))
User      System verstrichen
0.06        0.00        0.06
> all.equal(RES1, RES2)
[1] TRUE

There is about a 2-3 fold gain in speed for the above two setups but a 35-fold gain in the small window/large slides setup. Interesting…
I noticed that zoo:::rollapply.zoo uses mapply internally, maybe there is some overhead for pure vector calculations…

Cheers,
Andrej

## bigcor: Large correlation matrices in R

February 22, 2013

Here is the update to my bigcor function for creating large correlation/covariance matrices.
What does it do?
1) It splits a matrix with N columns into equal size blocks (+ a remainder block) of $A_1 ... A_n$$B_1 ... B_n$, etc. Block size can be defined by user, 2000 is a good value because cor can handle this quite quickly. If the matrix has 13796 columns, the split will be 2000;2000;2000;2000;2000;2000;1796.
2) For all combinations of blocks, the correlation matrix is calculated, so $latex A_n/A_n, A_n/B_n, B_n/B_n etc. 3) The blockwise correlation matrices are transferred into the preallocated empty matrix at the appropriate position (where the correlations would usually reside). To ensure symmetry around the diagonal, this is done twice in the upper and lower triangle. This way, the N x N empty correlation matrix is filled like a checkerboard with patches of n x n correlation sub-matrices. In our case, in which we split the N = 40000 columns into n = 8 blocks, we obtain 36 combinations (combn(1:8, 2) + 8; + 8 because of A/A, B/B etc) of blocks with dimension 5000 x 5000 each. This gives 36 x 5000 x 5000 x 2 (filling both triangles) – 8 x 5000 x 5000 (because the diagonals are created twice) = 1.6E9 = 40000 x 40000. The function can found here: http://www.dr-spiess.de/scripts/bigcor.r Timings for a 21796 x 21796 matrix roughly 2 min! ## The magic empty bracket January 30, 2013 I have been working with R for some time now, but once in a while, basic functions catch my eye that I was not aware of… For some project I wanted to transform a correlation matrix into a covariance matrix. Now, since cor2cov does not exist, I thought about “reversing” the cov2cor function (stats:::cov2cor). Inside the code of this function, a specific line jumped into my retina: r[] <- Is * V * rep(Is, each = p)  What’s this [ ]? Well, it stands for every element $E_{ij}$ of matrix $E$. Consider this: mat <- matrix(NA, nrow = 5, ncol = 5)  > mat [,1] [,2] [,3] [,4] [,5] [1,] NA NA NA NA NA [2,] NA NA NA NA NA [3,] NA NA NA NA NA [4,] NA NA NA NA NA [5,] NA NA NA NA NA With the empty bracket, we can now substitute ALL values by a new value: mat[] <- 1  > mat [,1] [,2] [,3] [,4] [,5] [1,] 1 1 1 1 1 [2,] 1 1 1 1 1 [3,] 1 1 1 1 1 [4,] 1 1 1 1 1 [5,] 1 1 1 1 1 Interestingly, this also works with lists: L <- list(a = 1, b = 2, c = 3) >L$a
[1] 1
$b [1] 2$c
[1] 3
L[] <- 5

> L
$a [1] 5$b
[1] 5

$c [1] 5 Cheers, Andrej  ## Peer-reviewed R packages? November 23, 2012 Dear R-Users, a question: I am the author of the ‘qpcR’ package. Within this, there is a function ‘propagate’ that does error propagation based on Monte Carlo Simulation, permutation-based confidence intervals and Taylor expansion. For the latter I recently implemented a second-order Taylor expansion term that can correct for nonlinearity. The formulas are quite complex and I had quite a hard time to transform the second-order term that has covariance included into a matrix form. The first order expansion can be represented by GRAD * SIGMA * t(GRAD) with GRAD being the gradient matrix and SIGMA the covariance matrix. The second order term is 0.5 * trace(HESS * SIGMA * HESS * SIGMA) with HESS being the Hessian matrix of partial second derivatives. So here the problem: Being a molecular biologist, my statistical skills have limitations. I am pretty sure my implementation is right (I checked the examples in GUM “Guide to the evaluation of uncertainties in experiments” and got the same results). BUT: Being sure says nothing. It could well be that I’m on the wrong track and people use a function that gives essentially wrong (maybe slightly wrong which may still be fatal…) results because they rely on me doing things right. Now let’s suppose I’m not the only one with a queer feeling when submitting packages to CRAN in the hope that everything I did is right. Which brings me to the question: What are the merits of a peer-review system for R packages? I would suppose (if the paradigm is right that peer-reviewing DOES increase quality…) that if packages were independently peer-reviewed in a way that a handful of other people inspect the code and check for detrimental bugs that are enigmatic to the packages author, then this would significantly increase the credibility of the packages functions. Now I know that there is the “R Journal” or the “Journal of Statistical Software” that regularly publish papers where new R packages are explained in detail. However, they are restricted to the statistical background of the implementations and I’m not sure if the lines of code are actually checked. An idea would be that such a review system is voluntary, but if a package has been reviewed and declared as “passed”, then this would be some sort of quality criterion. For scientists like me that need to write grant applications, this would give more impact on the work one has dedicated into developing a “good” package that could maybe be cited as being “peer-reviewed”. Pros: * increase of quality * maybe citable * beneficial for grant applications or C.V. Cons: * reviewers must meticulously inspect the code * longer time to CRAN submission * no platform yet to publish the reviews which can be cited like a journal Comments welcome! Cheers, Andrej ## A weighting function for ‘nls’ / ‘nlsLM’ July 19, 2012 Standard nonlinear regression assumes homoscedastic data, that is, all response values $y_i$ are distributed normally. In case of heteroscedastic data (i.e. when the variance is dependent on the magnitude of the data), weighting the fit is essential. In nls (or nlsLM of the minpack.lm package), weighting can be conducted by two different methods: 1) by supplying a vector of weighting values $w_i$ for every $y_i$ that is supplied to the ‘weights’ argument. 2) by modifying the objective function in a way that it has the weights included in the body of the function. An example is given for this one in the documentation to nls. 1) is easy of course, but all $w_i$ values have to be derived from somewhere, for example the reciproce of the response values $w_i = 1/y_i$. 2) is a bit more cumbersome and one does not always want to define a new weighted version of the objective function. I found it about time to design a weighting function that can be supplied to the ‘weights’ argument and that makes weighting by different regimes a breeze. wfct() returns a vector of weights that are calculated from a user-defined expression and transfers this vector within nls. wfct <- function(expr) { expr <- deparse(substitute(expr)) ## create new environment newEnv <- new.env() ## get call mc <- sys.calls()[[1]] mcL <- as.list(mc) ## get data and write to newEnv DATA <- mcL[["data"]] DATA <- eval(DATA) DATA <- as.list(DATA) NAMES <- names(DATA) for (i in 1:length(DATA)) assign(NAMES[i], DATA[[i]], envir = newEnv) ## get parameter, response and predictor names formula <- as.formula(mcL[[2]]) VARS <- all.vars(formula) RESP <- VARS[1] RHS <- VARS[-1] PRED <- match(RHS, names(DATA)) PRED <- names(DATA)[na.omit(PRED)] ## calculate variances for response values if "error" is in expression ## and write to newEnv if (length(grep("error", expr)) > 0) { y <- DATA[[RESP]] x <- DATA[[PRED]] ## test for replication if (!any(duplicated(x))) stop("No replicates available to calculate error from!") ## calculate error error <- tapply(y, x, function(e) var(e, na.rm = TRUE)) error <- as.numeric(sqrt(error)) ## convert to original repititions error <- rep(error, as.numeric(table(x))) assign("error", error, envir = newEnv) } ## calculate fitted or residual values if "fitted"/"resid" is in expression ## and write to newEnv if (length(grep("fitted", expr)) > 0 || length(grep("resid", expr)) > 0) { mc2 <- mc mc2$weights <- NULL
MODEL <- eval(mc2)
fitted <- fitted(MODEL)
resid <- residuals(MODEL)
assign("fitted", fitted, newEnv)
assign("resid", resid, newEnv)
}

## return evaluation in newEnv: vector of weights
OUT <- eval(parse(text = expr), envir = newEnv)
return(OUT)
}

The weighting function can take 5 different variable definitions and combinations thereof:
the name of the predictor (independent) variable,
the name of the response (dependent) variable,
error: if replicates $y_{ij}$ exist, the error $\sigma(y_{ij})$,
fitted: the fitted values $\hat{y}_i$ of the model,
resid: the residuals $y_i - \hat{y}_i$ of the model,
For the last two, the model is fit unweighted, fitted values and residuals are extracted and the model is refit by the defined weights.

Let’s have some examples (taken from the nls doc):

Treated <- Puromycin[Puromycin\$state == "treated", ]

Weighting by inverse of response $1/y_i$:

nls(rate ~ Vm * conc/(K + conc), data = Treated,
start = c(Vm = 200, K = 0.05), weights = wfct(1/rate))
Nonlinear regression model
model:  rate ~ Vm * conc/(K + conc)
data:  Treated
Vm         K
209.59669   0.06065
weighted residual sum-of-squares: 12.27

Number of iterations to convergence: 5
Achieved convergence tolerance: 5.2e-06

Weighting by square root of predictor $\sqrt{x_i}$:

nls(rate ~ Vm * conc/(K + conc), data = Treated,
start = c(Vm = 200, K = 0.05), weights = wfct(sqrt(conc)))
Nonlinear regression model
model:  rate ~ Vm * conc/(K + conc)
data:  Treated
Vm         K
216.74134   0.07106
weighted residual sum-of-squares: 332.6

Number of iterations to convergence: 5
Achieved convergence tolerance: 2.403e-06

Weighting by inverse square of fitted values $1/\hat{y_i}^2$:

nls(rate ~ Vm * conc/(K + conc), data = Treated,
start = c(Vm = 200, K = 0.05), weights = wfct(1/fitted^2))
Nonlinear regression model
model:  rate ~ Vm * conc/(K + conc)
data:  Treated
Vm         K
200.84098   0.04943
weighted residual sum-of-squares: 0.2272

Number of iterations to convergence: 4
Achieved convergence tolerance: 1.47e-06

A common weighting regime is to use the inverse variance $1/\sigma^2(y_i)$:

## A better ‘nls’ (?)

July 5, 2012

Those that do a lot of nonlinear regression will love the nls function of R. In most of the cases it works really well, but there are some mishaps that can occur when using bad starting values for the parameters. One of the most dreaded is the “singular gradient matrix at initial parameter estimates” which brings the function to a stop because the gradient check in stats:::nlsModel will terminate if the QR decomposition is not of full column rank.

Nearly all nonlinear fitting programs out there use the Levenberg-Marquardt algorithm for nonlinear regression. This is because its switching between Gauss-Newton and gradient descent is highly robust against far-off-optimal starting values.  Unfortunately, the standard nls function has no LM implemented, instead it houses the Gauss-Newton type, the PORT routines and a partial linear fitter. The fabulous minpack.lm package from Katherine M. Mullen offers an R frontend to a Fortran LM implementation of the MINPACK package. The nls.lm function must be supplied with an objective function that returns a vector of residuals to be minimized. Together with Kate I developed a function nlsLM that has the gee wizz formulaic interface of nls, but it calls LM instead of Gauss-Newton. This has some advantages as we will see below.  The function returns the usual result of class ‘nls’, and due to some modifications, all standard generics work. The modifications were made so that the formula is transformed into a function that returns a vector of (weighted) residuals whose sum square is minimized by nls.lm. The optimized parameters are then transferred to stats:::nlsModel in order to obtain an object of class ‘nlsModel’. The internal C  functions C_nls_iter and stats:::nls_port_fit were removed to avoid subsequent “Gauss-Newton”, “port” or “plinear” optimization of nlsModel.

So, what’s similar and what’s, well, better…

library(minpack.lm)
### Examples from 'nls' doc ###
DNase1 <- subset(DNase, Run == 1)
fm1DNase1 <- nlsLM(density ~ Asym/(1 + exp((xmid - log(conc))/scal)),
data = DNase1, start = list(Asym = 3, xmid = 0, scal = 1))
## all generics are applicable
coef(fm1DNase1)
confint(fm1DNase1)
deviance(fm1DNase1)
df.residual(fm1DNase1)
fitted(fm1DNase1)
formula(fm1DNase1)
logLik(fm1DNase1)
predict(fm1DNase1)
print(fm1DNase1)
profile(fm1DNase1)
residuals(fm1DNase1)
summary(fm1DNase1)
update(fm1DNase1)
vcov(fm1DNase1)
weights(fm1DNase1)

nlsLM can fit zero noise data, when nls fails:
x <- 1:10
y <- 2*x + 3
nls(y ~ a + b * x, start = list(a = 0.12345, b = 0.54321))
nlsLM(y ~ a + b * x, start = list(a = 0.12345, b = 0.54321))

Nonlinear regression model
model:  y ~ a + b * x
data:  parent.frame()
a b
3 2
residual sum-of-squares: 0
Number of iterations to convergence: 3
Achieved convergence tolerance: 1.49e-08
Example taken from here:
x <- 0:140
y <- 200 / (1 + exp(17 - x)/2) * exp(-0.02*x)
yeps <- y + rnorm(length(y), sd = 2)
nls(yeps ~ p1 / (1 + exp(p2 - x)) * exp(p4 * x), start=list(p1=410,p2=18,p4=-.03))

Optimal starting parameters work:
Nonlinear regression model
model:  yeps ~ p1/(1 + exp(p2 - x)) * exp(p4 * x)
data:  parent.frame()
p1        p2        p4
200.67285  16.32430  -0.02005
residual sum-of-squares: 569
Number of iterations to convergence: 6
Achieved convergence tolerance: 2.611e-07
But off-optimal parameters give error:
nls(yeps ~ p1 / (1 + exp(p2 - x)) * exp(p4 * x), start=list(p1=10,p2=18,p4=-.03))

Fehler in nls(yeps ~ p1/(1 + exp(p2 - x)) * exp(p4 * x), start = list(p1 = 10,  :
nlsLM is more robust with these starting parameters:
nlsLM(yeps ~ p1 / (1 + exp(p2 - x)) * exp(p4 * x), start=list(p1=10,p2=18,p4=-.03))

Nonlinear regression model
model:  yeps ~ p1/(1 + exp(p2 - x)) * exp(p4 * x)
data:  parent.frame()
p1        p2        p4
200.67285  16.32430  -0.02005
residual sum-of-squares: 569
Number of iterations to convergence: 10
Achieved convergence tolerance: 1.49e-08
Have fun!
Cheers, Andrej


text:
If the sum of three consecutive integers is 9, what are the integers?

Oct 26, 2016

Consecutive integers are $2$, $3$ and $4$

Explanation:

Let the consecutive integers be $x$, $x + 1$ and $x + 2$.

As their sum is $9$, we have

$x + x + 1 + x + 2 = 9$

or $3 x + 3 = 9$

or $3 x = 9 - 3 = 6$

Hence, x=6×1/3=2

and consecutive integers are $2$, $3$ and $4$.