Symbolic Regression of Integer Sequences (SRIS) aims to discover precise mathematical formulas from integer sequences. The neural machine translation-based method of SRIS trains the model using randomly generated data, and directly utilizes the trained model for inference on target sequences. However, the method often fails to effectively generalize to the target sequence, since the randomly generated data can not adequately cover the distributions of target data, i.e., there are distribution differences between them. In this work, we propose a progressive self-learning (PSL) method to explicitly capture sequence-formula distributions of the target domain. Specifically, a source domain dataset is generated by incorporating initial terms of the target domain to reduce the sequence distribution gap between the source domain and the target domain. Meanwhile, a self-learning loop strategy is adopted to improve the ability of the model to capture the sequence-formula distribution of the target domain. In this strategy, a neural machine translation model is used to learn the mappings from sequences to formulas in an end-to-end fashion. Then, this model is employed to explore candidate formulas of the target sequence using beam search. After verifying these candidate formula correctness, some of them are retained as training data for the next learning. Experimental results on OEIS datasets demonstrate that the proposed method surpasses current state-of-the-art methods in accuracy, and also discovers new formulas.